MFEM  v3.3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
triangle.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "mesh_headers.hpp"
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 
37 int Triangle::NeedRefinement(DSTable &v_to_v, int *middle) const
38 {
39  int m;
40 
41  if ((m = v_to_v(indices[0], indices[1])) != -1 && middle[m] != -1) { return 1; }
42  if ((m = v_to_v(indices[1], indices[2])) != -1 && middle[m] != -1) { return 1; }
43  if ((m = v_to_v(indices[2], indices[0])) != -1 && middle[m] != -1) { return 1; }
44  return 0;
45 }
46 
47 void Triangle::SetVertices(const int *ind)
48 {
49  for (int i = 0; i < 3; i++)
50  {
51  indices[i] = ind[i];
52  }
53 }
54 
56 {
57  double d[3];
58  int shift, v;
59 
60  d[0] = ( (pmat(0,1)-pmat(0,0))*(pmat(0,1)-pmat(0,0)) +
61  (pmat(1,1)-pmat(1,0))*(pmat(1,1)-pmat(1,0)) );
62  d[1] = ( (pmat(0,2)-pmat(0,1))*(pmat(0,2)-pmat(0,1)) +
63  (pmat(1,2)-pmat(1,1))*(pmat(1,2)-pmat(1,1)) );
64  d[2] = ( (pmat(0,2)-pmat(0,0))*(pmat(0,2)-pmat(0,0)) +
65  (pmat(1,2)-pmat(1,0))*(pmat(1,2)-pmat(1,0)) );
66 
67  // if pmat has 3 rows, then use extra term in each sum
68  if (pmat.Height()==3)
69  {
70  d[0] += (pmat(2,1)-pmat(2,0))*(pmat(2,1)-pmat(2,0));
71  d[1] += (pmat(2,2)-pmat(2,1))*(pmat(2,2)-pmat(2,1));
72  d[2] += (pmat(2,2)-pmat(2,0))*(pmat(2,2)-pmat(2,0));
73  }
74 
75  if (d[0] >= d[1])
76  {
77  if (d[0] >= d[2]) { shift = 0; }
78  else { shift = 2; }
79  }
80  else if (d[1] >= d[2]) { shift = 1; }
81  else { shift = 2; }
82 
83  switch (shift)
84  {
85  case 0:
86  break;
87  case 1:
88  v = indices[0];
89  indices[0] = indices[1];
90  indices[1] = indices[2];
91  indices[2] = v;
92  break;
93  case 2:
94  v = indices[0];
95  indices[0] = indices[2];
96  indices[2] = indices[1];
97  indices[1] = v;
98  break;
99  }
100 }
101 
102 void Triangle::MarkEdge(const DSTable &v_to_v, const int *length)
103 {
104  int l, L, j, ind[3], i;
105 
106  L = length[ v_to_v(indices[0], indices[1]) ]; j = 0;
107  if ( (l = length[ v_to_v(indices[1], indices[2]) ]) > L ) { L = l; j = 1; }
108  if ( (l = length[ v_to_v(indices[2], indices[0]) ]) > L ) { L = l; j = 2; }
109 
110  for (i = 0; i < 3; i++)
111  {
112  ind[i] = indices[i];
113  }
114 
115  switch (j)
116  {
117  case 1:
118  indices[0] = ind[1]; indices[1] = ind[2]; indices[2] = ind[0];
119  break;
120  case 2:
121  indices[0] = ind[2]; indices[1] = ind[0]; indices[2] = ind[1];
122  break;
123  }
124 }
125 
126 // static method
127 void Triangle::GetPointMatrix(unsigned transform, DenseMatrix &pm)
128 {
129  double *a = &pm(0,0), *b = &pm(0,1), *c = &pm(0,2);
130 
131  // initialize to identity
132  a[0] = 0.0; a[1] = 0.0;
133  b[0] = 1.0; b[1] = 0.0;
134  c[0] = 0.0; c[1] = 1.0;
135 
136  int chain[12], n = 0;
137  while (transform)
138  {
139  chain[n++] = (transform & 7) - 1;
140  transform >>= 3;
141  }
142 
143  /* The transformations and orientations here match
144  Mesh::UniformRefinement and Mesh::Bisection for triangles:
145 
146  c c
147  * *
148  | \ |\\
149  | \ | \ \
150  | 2 \ e | \ \
151  f *-------* | \ \
152  | \ 3 | \ | \ \
153  | \ | \ | 4 \ 5 \
154  | 0 \ | 1 \ | \ \
155  *-------*-------* *-------*-------*
156  a d b a d b
157  */
158 
159  double d[2], e[2], f[2];
160 #define ASGN(a, b) (a[0] = b[0], a[1] = b[1])
161 #define AVG(a, b, c) (a[0] = (b[0] + c[0])*0.5, a[1] = (b[1] + c[1])*0.5)
162 
163  while (n)
164  {
165  switch (chain[--n])
166  {
167  case 0: AVG(b, a, b); AVG(c, a, c); break;
168  case 1: AVG(a, a, b); AVG(c, b, c); break;
169  case 2: AVG(a, a, c); AVG(b, b, c); break;
170 
171  case 3:
172  AVG(d, a, b); AVG(e, b, c); AVG(f, c, a);
173  ASGN(a, e); ASGN(b, f); ASGN(c, d); break;
174 
175  case 4:
176  AVG(d, a, b); // N.B.: orientation
177  ASGN(b, a); ASGN(a, c); ASGN(c, d); break;
178 
179  case 5:
180  AVG(d, a, b); // N.B.: orientation
181  ASGN(a, b); ASGN(b, c); ASGN(c, d); break;
182 
183  default:
184  MFEM_ABORT("Invalid transform.");
185  }
186  }
187 }
188 
190 {
191  v.SetSize(3);
192  for (int i = 0; i < 3; i++)
193  {
194  v[i] = indices[i];
195  }
196 }
197 
199 
200 } // namespace mfem
virtual int * GetVertices()
Definition: triangle.hpp:71
virtual int NeedRefinement(DSTable &v_to_v, int *middle) const
Return 1 if the element needs refinement in order to get conforming mesh.
Definition: triangle.cpp:37
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void MarkEdge(DenseMatrix &pmat)
Definition: triangle.cpp:55
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
unsigned transform
Definition: triangle.hpp:28
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
Definition: triangle.cpp:47
Class for linear FE on triangle.
Definition: fe.hpp:520
int attribute
Element&#39;s attribute (specifying material property, etc).
Definition: element.hpp:32
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:503
Linear2DFiniteElement TriangleFE
Definition: triangle.cpp:198
Abstract data type element.
Definition: element.hpp:27
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
Definition: triangle.cpp:127
int indices[3]
Definition: triangle.hpp:26