MFEM  v4.1.0 Finite element discretization library
triangle.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
13
14 namespace mfem
15 {
16
17 Triangle::Triangle(const int *ind, int attr) : Element(Geometry::TRIANGLE)
18 {
19  attribute = attr;
20  for (int i = 0; i < 3; i++)
21  {
22  indices[i] = ind[i];
23  }
24  transform = 0;
25 }
26
27 Triangle::Triangle(int ind1, int ind2, int ind3, int attr)
28  : Element(Geometry::TRIANGLE)
29 {
30  attribute = attr;
31  indices[0] = ind1;
32  indices[1] = ind2;
33  indices[2] = ind3;
34  transform = 0;
35 }
36
38 {
39  if (v_to_v.FindId(indices[0], indices[1]) != -1) { return 1; }
40  if (v_to_v.FindId(indices[1], indices[2]) != -1) { return 1; }
41  if (v_to_v.FindId(indices[2], indices[0]) != -1) { return 1; }
42  return 0;
43 }
44
45 void Triangle::SetVertices(const int *ind)
46 {
47  for (int i = 0; i < 3; i++)
48  {
49  indices[i] = ind[i];
50  }
51 }
52
54 {
55  double d[3];
56  int shift, v;
57
58  d[0] = ( (pmat(0,1)-pmat(0,0))*(pmat(0,1)-pmat(0,0)) +
59  (pmat(1,1)-pmat(1,0))*(pmat(1,1)-pmat(1,0)) );
60  d[1] = ( (pmat(0,2)-pmat(0,1))*(pmat(0,2)-pmat(0,1)) +
61  (pmat(1,2)-pmat(1,1))*(pmat(1,2)-pmat(1,1)) );
62  d[2] = ( (pmat(0,2)-pmat(0,0))*(pmat(0,2)-pmat(0,0)) +
63  (pmat(1,2)-pmat(1,0))*(pmat(1,2)-pmat(1,0)) );
64
65  // if pmat has 3 rows, then use extra term in each sum
66  if (pmat.Height()==3)
67  {
68  d[0] += (pmat(2,1)-pmat(2,0))*(pmat(2,1)-pmat(2,0));
69  d[1] += (pmat(2,2)-pmat(2,1))*(pmat(2,2)-pmat(2,1));
70  d[2] += (pmat(2,2)-pmat(2,0))*(pmat(2,2)-pmat(2,0));
71  }
72
73  if (d[0] >= d[1])
74  {
75  if (d[0] >= d[2]) { shift = 0; }
76  else { shift = 2; }
77  }
78  else if (d[1] >= d[2]) { shift = 1; }
79  else { shift = 2; }
80
81  switch (shift)
82  {
83  case 0:
84  break;
85  case 1:
86  v = indices[0];
87  indices[0] = indices[1];
88  indices[1] = indices[2];
89  indices[2] = v;
90  break;
91  case 2:
92  v = indices[0];
93  indices[0] = indices[2];
94  indices[2] = indices[1];
95  indices[1] = v;
96  break;
97  }
98 }
99
100 // Static method
101 void Triangle::MarkEdge(int *indices, const DSTable &v_to_v, const int *length)
102 {
103  int l, L, j, ind[3], i;
104
105  L = length[ v_to_v(indices[0], indices[1]) ]; j = 0;
106  if ( (l = length[ v_to_v(indices[1], indices[2]) ]) > L ) { L = l; j = 1; }
107  if ( (l = length[ v_to_v(indices[2], indices[0]) ]) > L ) { j = 2; }
108
109  for (i = 0; i < 3; i++)
110  {
111  ind[i] = indices[i];
112  }
113
114  switch (j)
115  {
116  case 1:
117  indices[0] = ind[1]; indices[1] = ind[2]; indices[2] = ind[0];
118  break;
119  case 2:
120  indices[0] = ind[2]; indices[1] = ind[0]; indices[2] = ind[1];
121  break;
122  }
123 }
124
125 // static method
126 void Triangle::GetPointMatrix(unsigned transform, DenseMatrix &pm)
127 {
128  double *a = &pm(0,0), *b = &pm(0,1), *c = &pm(0,2);
129
130  // initialize to identity
131  a[0] = 0.0; a[1] = 0.0;
132  b[0] = 1.0; b[1] = 0.0;
133  c[0] = 0.0; c[1] = 1.0;
134
135  int chain[12], n = 0;
136  while (transform)
137  {
138  chain[n++] = (transform & 7) - 1;
139  transform >>= 3;
140  }
141
142  /* The transformations and orientations here match
143  Mesh::UniformRefinement and Mesh::Bisection for triangles:
144
145  c c
146  * *
147  | \ |\\
148  | \ | \ \
149  | 2 \ e | \ \
150  f *-------* | \ \
151  | \ 3 | \ | \ \
152  | \ | \ | 4 \ 5 \
153  | 0 \ | 1 \ | \ \
154  *-------*-------* *-------*-------*
155  a d b a d b
156  */
157
158  double d[2], e[2], f[2];
159 #define ASGN(a, b) (a[0] = b[0], a[1] = b[1])
160 #define AVG(a, b, c) (a[0] = (b[0] + c[0])*0.5, a[1] = (b[1] + c[1])*0.5)
161
162  while (n)
163  {
164  switch (chain[--n])
165  {
166  case 0: AVG(b, a, b); AVG(c, a, c); break;
167  case 1: AVG(a, a, b); AVG(c, b, c); break;
168  case 2: AVG(a, a, c); AVG(b, b, c); break;
169
170  case 3:
171  AVG(d, a, b); AVG(e, b, c); AVG(f, c, a);
172  ASGN(a, e); ASGN(b, f); ASGN(c, d); break;
173
174  case 4:
175  AVG(d, a, b); // N.B.: orientation
176  ASGN(b, a); ASGN(a, c); ASGN(c, d); break;
177
178  case 5:
179  AVG(d, a, b); // N.B.: orientation
180  ASGN(a, b); ASGN(b, c); ASGN(c, d); break;
181
182  default:
183  MFEM_ABORT("Invalid transform.");
184  }
185  }
186 }
187
189 {
190  v.SetSize(3);
191  for (int i = 0; i < 3; i++)
192  {
193  v[i] = indices[i];
194  }
195 }
196
197 } // namespace mfem
virtual int * GetVertices()
Definition: triangle.hpp:74
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void MarkEdge(DenseMatrix &pmat)
Definition: triangle.cpp:53
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
double b
Definition: lissajous.cpp:42
unsigned transform
Definition: triangle.hpp:28
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
Definition: triangle.cpp:45
int attribute
Element&#39;s attribute (specifying material property, etc).
Definition: element.hpp:33
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
double a
Definition: lissajous.cpp:41
virtual int NeedRefinement(HashTable< Hashed2 > &v_to_v) const
Return 1 if the element needs refinement in order to get conforming mesh.
Definition: triangle.cpp:37
int FindId(int p1, int p2) const
Find id of item whose parents are p1, p2... Return -1 if it doesn&#39;t exist.
Definition: hash.hpp:370
Abstract data type element.
Definition: element.hpp:28
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
Definition: triangle.cpp:126
int indices[3]
Definition: triangle.hpp:26