gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
MTetrahedron.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include "GmshConfig.h"
7 #include "MTetrahedron.h"
8 #include "Numeric.h"
9 #include "Context.h"
10 #include "BasisFactory.h"
11 #include "pointsGenerators.h"
12 
13 #if defined(HAVE_MESH)
14 #include "qualityMeasures.h"
15 #include "meshGFaceDelaunayInsertion.h"
16 #include "meshGRegionDelaunayInsertion.h"
17 #endif
18 
19 #define SQU(a) ((a) * (a))
20 
21 std::map<int, IndicesReversed> MTetrahedronN::_order2indicesReversedTet;
22 
23 void MTetrahedron::getEdgeRep(bool curved, int num, double *x, double *y,
24  double *z, SVector3 *n)
25 {
26  // don't use MElement::_getEdgeRep: it's slow due to the creation of MFaces
27  MVertex *v0 = _v[edges_tetra(num, 0)];
28  MVertex *v1 = _v[edges_tetra(num, 1)];
29  x[0] = v0->x();
30  y[0] = v0->y();
31  z[0] = v0->z();
32  x[1] = v1->x();
33  y[1] = v1->y();
34  z[1] = v1->z();
35  if(CTX::instance()->mesh.lightLines > 1) {
36  static const int vv[6] = {2, 0, 1, 1, 0, 2};
37  MVertex *v2 = _v[vv[num]];
38  SVector3 t1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
39  SVector3 t2(v2->x() - x[0], v2->y() - y[0], v2->z() - z[0]);
40  SVector3 normal = crossprod(t1, t2);
41  normal.normalize();
42  n[0] = n[1] = normal;
43  }
44  else {
45  n[0] = n[1] = SVector3(0., 0., 1.);
46  }
47 }
48 
49 SPoint3 MTetrahedron::circumcenter()
50 {
51 #if defined(HAVE_MESH)
52  MTet4 t(this, 0);
53  double res[3];
54  t.circumcenter(res);
55  return SPoint3(res[0], res[1], res[2]);
56 #else
57  return SPoint3(0., 0., 0.);
58 #endif
59 }
60 
61 double MTetrahedron::getOuterRadius()
62 {
63 #if defined(HAVE_MESH)
64  SPoint3 center = circumcenter();
65  const double dx = getVertex(0)->x() - center.x();
66  const double dy = getVertex(0)->y() - center.y();
67  const double dz = getVertex(0)->z() - center.z();
68  double circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
69  return circum_radius;
70 #else
71  return 0.0;
72 #endif
73 }
74 
75 double MTetrahedron::getInnerRadius()
76 {
77  // radius of inscribed sphere = 3 * Volume / sum(Area_i)
78  double dist[3], face_area = 0.;
79  double vol = getVolume();
80  for(int i = 0; i < 4; i++) {
81  MFace f = getFace(i);
82  for(int j = 0; j < 3; j++) {
83  MEdge e = f.getEdge(j);
84  dist[j] = e.getVertex(0)->distance(e.getVertex(1));
85  }
86  face_area +=
87  0.25 *
88  sqrt((dist[0] + dist[1] + dist[2]) * (-dist[0] + dist[1] + dist[2]) *
89  (dist[0] - dist[1] + dist[2]) * (dist[0] + dist[1] - dist[2]));
90  }
91  return 3 * vol / face_area;
92 }
93 
94 double MTetrahedron::gammaShapeMeasure()
95 {
96 #if defined(HAVE_MESH)
97  double vol;
98  return qmTetrahedron::qm(this, qmTetrahedron::QMTET_GAMMA, &vol);
99 #else
100  return 0.;
101 #endif
102 }
103 
104 double MTetrahedron::etaShapeMeasure()
105 {
106 #if defined(HAVE_MESH)
107  double vol;
108  return qmTetrahedron::qm(this, qmTetrahedron::QMTET_ETA, &vol);
109 #else
110  return 0.;
111 #endif
112 }
113 
114 double MTetrahedron::getVolume()
115 {
116  double mat[3][3];
117  getMat(mat);
118  return det3x3(mat) / 6.;
119 }
120 
121 void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3]) const
122 {
123  double mat[3][3], b[3], det;
124  getMat(mat);
125  b[0] = xyz[0] - getVertex(0)->x();
126  b[1] = xyz[1] - getVertex(0)->y();
127  b[2] = xyz[2] - getVertex(0)->z();
128  sys3x3(mat, b, uvw, &det);
129 }
130 
131 int MTetrahedron::numCommonNodesInDualGraph(const MElement *const other) const
132 {
133  switch(other->getType()) {
134  case TYPE_PNT: return 1;
135  case TYPE_LIN: return 2;
136  case TYPE_QUA: return 4;
137  default: return 3;
138  }
139 }
140 
141 int MTetrahedron10::getNumEdgesRep(bool curved)
142 {
143  return curved ? 6 * CTX::instance()->mesh.numSubEdges : 6;
144 }
145 
146 int MTetrahedronN::getNumEdgesRep(bool curved)
147 {
148  return curved ? 6 * CTX::instance()->mesh.numSubEdges : 6;
149 }
150 
151 static void _myGetEdgeRep(MTetrahedron *tet, int num, double *x, double *y,
152  double *z, SVector3 *n, int numSubEdges)
153 {
154  static double pp[4][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
155  static int ed[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}};
156  int iEdge = num / numSubEdges;
157  int iSubEdge = num % numSubEdges;
158 
159  int iVertex1 = ed[iEdge][0];
160  int iVertex2 = ed[iEdge][1];
161  double t1 = (double)iSubEdge / (double)numSubEdges;
162  double u1 = pp[iVertex1][0] * (1. - t1) + pp[iVertex2][0] * t1;
163  double v1 = pp[iVertex1][1] * (1. - t1) + pp[iVertex2][1] * t1;
164  double w1 = pp[iVertex1][2] * (1. - t1) + pp[iVertex2][2] * t1;
165 
166  double t2 = (double)(iSubEdge + 1) / (double)numSubEdges;
167  double u2 = pp[iVertex1][0] * (1. - t2) + pp[iVertex2][0] * t2;
168  double v2 = pp[iVertex1][1] * (1. - t2) + pp[iVertex2][1] * t2;
169  double w2 = pp[iVertex1][2] * (1. - t2) + pp[iVertex2][2] * t2;
170 
171  SPoint3 pnt1, pnt2;
172  tet->pnt(u1, v1, w1, pnt1);
173  tet->pnt(u2, v2, w2, pnt2);
174  x[0] = pnt1.x();
175  x[1] = pnt2.x();
176  y[0] = pnt1.y();
177  y[1] = pnt2.y();
178  z[0] = pnt1.z();
179  z[1] = pnt2.z();
180 
181  // not great, but better than nothing
182  static const int f[6] = {0, 0, 0, 1, 2, 3};
183  n[0] = n[1] = tet->getFace(f[iEdge]).normal();
184 }
185 
186 void MTetrahedron10::getEdgeRep(bool curved, int num, double *x, double *y,
187  double *z, SVector3 *n)
188 {
189  if(curved)
190  _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
191  else
192  MTetrahedron::getEdgeRep(false, num, x, y, z, n);
193 }
194 
195 void MTetrahedronN::getEdgeRep(bool curved, int num, double *x, double *y,
196  double *z, SVector3 *n)
197 {
198  if(curved)
199  _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
200  else
201  MTetrahedron::getEdgeRep(false, num, x, y, z, n);
202 }
203 
204 int MTetrahedronN::getNumFacesRep(bool curved)
205 {
206  return curved ? 4 * SQU(CTX::instance()->mesh.numSubEdges) : 4;
207 }
208 
209 int MTetrahedron10::getNumFacesRep(bool curved)
210 {
211  return curved ? 4 * SQU(CTX::instance()->mesh.numSubEdges) : 4;
212 }
213 
214 static void _myGetFaceRep(MTetrahedron *tet, int num, double *x, double *y,
215  double *z, SVector3 *n, int numSubEdges)
216 {
217  static double pp[4][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
218  static int fak[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
219  int iFace = num / (numSubEdges * numSubEdges);
220  int iSubFace = num % (numSubEdges * numSubEdges);
221 
222  int iVertex1 = fak[iFace][0];
223  int iVertex2 = fak[iFace][1];
224  int iVertex3 = fak[iFace][2];
225 
226  /*
227  0
228  0 1
229  0 1 2
230  0 1 2 3
231  0 1 2 3 4
232  0 1 2 3 4 5
233  */
234 
235  // on the first layer, we have (numSubEdges-1) * 2 + 1 triangles
236  // on the second layer, we have (numSubEdges-2) * 2 + 1 triangles
237  // on the ith layer, we have (numSubEdges-1-i) * 2 + 1 triangles
238  int ix = 0, iy = 0;
239  int nbt = 0;
240  for(int i = 0; i < numSubEdges; i++) {
241  int nbl = (numSubEdges - i - 1) * 2 + 1;
242  nbt += nbl;
243  if(nbt > iSubFace) {
244  iy = i;
245  ix = nbl - (nbt - iSubFace);
246  break;
247  }
248  }
249 
250  const double d = 1. / numSubEdges;
251 
252  SPoint3 pnt1, pnt2, pnt3;
253  double u1, v1, u2, v2, u3, v3;
254  if(ix % 2 == 0) {
255  u1 = ix / 2 * d;
256  v1 = iy * d;
257  u2 = (ix / 2 + 1) * d;
258  v2 = iy * d;
259  u3 = ix / 2 * d;
260  v3 = (iy + 1) * d;
261  }
262  else {
263  u1 = (ix / 2 + 1) * d;
264  v1 = iy * d;
265  u2 = (ix / 2 + 1) * d;
266  v2 = (iy + 1) * d;
267  u3 = ix / 2 * d;
268  v3 = (iy + 1) * d;
269  }
270 
271  double U1 = pp[iVertex1][0] * (1. - u1 - v1) + pp[iVertex2][0] * u1 +
272  pp[iVertex3][0] * v1;
273  double U2 = pp[iVertex1][0] * (1. - u2 - v2) + pp[iVertex2][0] * u2 +
274  pp[iVertex3][0] * v2;
275  double U3 = pp[iVertex1][0] * (1. - u3 - v3) + pp[iVertex2][0] * u3 +
276  pp[iVertex3][0] * v3;
277 
278  double V1 = pp[iVertex1][1] * (1. - u1 - v1) + pp[iVertex2][1] * u1 +
279  pp[iVertex3][1] * v1;
280  double V2 = pp[iVertex1][1] * (1. - u2 - v2) + pp[iVertex2][1] * u2 +
281  pp[iVertex3][1] * v2;
282  double V3 = pp[iVertex1][1] * (1. - u3 - v3) + pp[iVertex2][1] * u3 +
283  pp[iVertex3][1] * v3;
284 
285  double W1 = pp[iVertex1][2] * (1. - u1 - v1) + pp[iVertex2][2] * u1 +
286  pp[iVertex3][2] * v1;
287  double W2 = pp[iVertex1][2] * (1. - u2 - v2) + pp[iVertex2][2] * u2 +
288  pp[iVertex3][2] * v2;
289  double W3 = pp[iVertex1][2] * (1. - u3 - v3) + pp[iVertex2][2] * u3 +
290  pp[iVertex3][2] * v3;
291 
292  tet->pnt(U1, V1, W1, pnt1);
293  tet->pnt(U2, V2, W2, pnt2);
294  tet->pnt(U3, V3, W3, pnt3);
295 
296  x[0] = pnt1.x();
297  x[1] = pnt2.x();
298  x[2] = pnt3.x();
299  y[0] = pnt1.y();
300  y[1] = pnt2.y();
301  y[2] = pnt3.y();
302  z[0] = pnt1.z();
303  z[1] = pnt2.z();
304  z[2] = pnt3.z();
305 
306  SVector3 d1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
307  SVector3 d2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
308  n[0] = crossprod(d1, d2);
309  n[0].normalize();
310  n[1] = n[0];
311  n[2] = n[0];
312 }
313 
314 void MTetrahedronN::getFaceRep(bool curved, int num, double *x, double *y,
315  double *z, SVector3 *n)
316 {
317  if(curved)
318  _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
319  else
320  MTetrahedron::getFaceRep(false, num, x, y, z, n);
321 }
322 
323 void MTetrahedron10::getFaceRep(bool curved, int num, double *x, double *y,
324  double *z, SVector3 *n)
325 {
326  if(curved)
327  _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
328  else
329  MTetrahedron::getFaceRep(false, num, x, y, z, n);
330 }
331 
332 void MTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
333 {
334  *npts = getNGQTetPts(pOrder);
335  *pts = getGQTetPts(pOrder);
336 }
337 
338 bool MTetrahedron::getFaceInfo(const MFace &face, int &ithFace, int &sign,
339  int &rot) const
340 {
341  for(ithFace = 0; ithFace < 4; ithFace++) {
342  if(_getFaceInfo(getFace(ithFace), face, sign, rot)) return true;
343  }
344  Msg::Error("Could not get face information for tetrahedron %d", getNum());
345  return false;
346 }
347 
348 void _getIndicesReversedTet(int order, IndicesReversed &indices)
349 {
350  fullMatrix<double> ref = gmshGenerateMonomialsTetrahedron(order);
351 
352  indices.resize(ref.size1());
353  for(int i = 0; i < ref.size1(); ++i) {
354  const double u = ref(i, 0);
355  const double v = ref(i, 1);
356  const double w = ref(i, 2);
357  for(int j = 0; j < ref.size1(); ++j) {
358  if(u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
359  indices[i] = j;
360  break;
361  }
362  }
363  }
364 }
365 
366 void MTetrahedronN::reverse()
367 {
368  auto it = _order2indicesReversedTet.find(_order);
369  if(it == _order2indicesReversedTet.end()) {
370  IndicesReversed indices;
371  _getIndicesReversedTet(_order, indices);
372  _order2indicesReversedTet[_order] = indices;
373  it = _order2indicesReversedTet.find(_order);
374  }
375 
376  IndicesReversed &indices = it->second;
377 
378  // copy vertices
379  std::vector<MVertex *> oldv(4 + _vs.size());
380  std::copy(_v, _v + 4, oldv.begin());
381  std::copy(_vs.begin(), _vs.end(), oldv.begin() + 4);
382 
383  // reverse
384  for(int i = 0; i < 4; ++i) { _v[i] = oldv[indices[i]]; }
385  for(std::size_t i = 0; i < _vs.size(); ++i) { _vs[i] = oldv[indices[4 + i]]; }
386 }
MElement::getNum
virtual std::size_t getNum() const
Definition: MElement.h:68
crossprod
SVector3 crossprod(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:150
MTetrahedronN::_order2indicesReversedTet
static std::map< int, IndicesReversed > _order2indicesReversedTet
Definition: MTetrahedron.h:373
MTetrahedronN::getFaceRep
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MTetrahedron.cpp:314
qualityMeasures.h
MTetrahedron::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3]) const
Definition: MTetrahedron.cpp:121
MEdge
Definition: MEdge.h:14
meshGRegionDelaunayInsertion.h
MTetrahedron::getOuterRadius
virtual double getOuterRadius()
Definition: MTetrahedron.cpp:61
TYPE_LIN
#define TYPE_LIN
Definition: GmshDefines.h:65
MTetrahedron
Definition: MTetrahedron.h:34
MTetrahedronN::_order
const char _order
Definition: MTetrahedron.h:377
MTetrahedron10::getNumEdgesRep
virtual int getNumEdgesRep(bool curved)
Definition: MTetrahedron.cpp:141
MTetrahedron::circumcenter
virtual SPoint3 circumcenter()
Definition: MTetrahedron.cpp:49
MVertex
Definition: MVertex.h:24
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
det3x3
double det3x3(double mat[3][3])
Definition: Numeric.cpp:126
SPoint3
Definition: SPoint3.h:14
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
MFace::normal
SVector3 normal() const
Definition: MFace.cpp:204
TYPE_PNT
#define TYPE_PNT
Definition: GmshDefines.h:64
SVector3
Definition: SVector3.h:16
getGQTetPts
IntPt * getGQTetPts(int order, bool forceTensorRule=false)
Definition: GaussQuadratureTet.cpp:3346
MTetrahedron::_v
MVertex * _v[4]
Definition: MTetrahedron.h:36
MTetrahedron::getInnerRadius
virtual double getInnerRadius()
Definition: MTetrahedron.cpp:75
w1
const double w1
Definition: GaussQuadratureHex.cpp:18
sys3x3
int sys3x3(double mat[3][3], double b[3], double res[3], double *det)
Definition: Numeric.cpp:147
MTetrahedronN::_vs
std::vector< MVertex * > _vs
Definition: MTetrahedron.h:376
MTet4
Definition: meshGRegionDelaunayInsertion.h:46
MTetrahedron::getMat
void getMat(double mat[3][3]) const
Definition: MTetrahedron.h:124
MTetrahedron::getEdgeRep
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MTetrahedron.cpp:23
CTX::instance
static CTX * instance()
Definition: Context.cpp:122
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
fullMatrix< double >
_myGetEdgeRep
static void _myGetEdgeRep(MTetrahedron *tet, int num, double *x, double *y, double *z, SVector3 *n, int numSubEdges)
Definition: MTetrahedron.cpp:151
MFace
Definition: MFace.h:20
MEdge::getVertex
MVertex * getVertex(std::size_t i) const
Definition: MEdge.h:39
contextMeshOptions::numSubEdges
int numSubEdges
Definition: Context.h:85
MElement::getType
virtual int getType() const =0
MTetrahedron::numCommonNodesInDualGraph
virtual int numCommonNodesInDualGraph(const MElement *const other) const
Definition: MTetrahedron.cpp:131
_getIndicesReversedTet
void _getIndicesReversedTet(int order, IndicesReversed &indices)
Definition: MTetrahedron.cpp:348
MTetrahedron::getVertex
virtual MVertex * getVertex(int num)
Definition: MTetrahedron.h:67
getNGQTetPts
int getNGQTetPts(int order, bool forceTensorRule=false)
Definition: GaussQuadratureTet.cpp:3362
Numeric.h
MTetrahedron::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MTetrahedron.cpp:332
MTetrahedron::getFaceInfo
virtual bool getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) const
Definition: MTetrahedron.cpp:338
MTetrahedronN::getNumFacesRep
virtual int getNumFacesRep(bool curved)
Definition: MTetrahedron.cpp:204
qmTetrahedron::qm
static double qm(MTetrahedron *t, const Measures &cr, double *volume=nullptr)
Definition: qualityMeasures.cpp:616
meshGFaceDelaunayInsertion.h
MElement::_getFaceInfo
static bool _getFaceInfo(const MFace &face, const MFace &other, int &sign, int &rot)
Definition: MElement.cpp:66
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
qmTetrahedron::QMTET_ETA
@ QMTET_ETA
Definition: qualityMeasures.h:68
CTX::mesh
contextMeshOptions mesh
Definition: Context.h:313
MElement
Definition: MElement.h:30
MTetrahedron::getVolume
virtual double getVolume()
Definition: MTetrahedron.cpp:114
MTetrahedron::edges_tetra
static int edges_tetra(const int edge, const int vert)
Definition: MTetrahedron.h:183
MElement::pnt
virtual void pnt(double u, double v, double w, SPoint3 &p) const
Definition: MElement.cpp:1072
TYPE_QUA
#define TYPE_QUA
Definition: GmshDefines.h:67
MTetrahedron10::getFaceRep
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MTetrahedron.cpp:323
Context.h
IntPt
Definition: GaussIntegration.h:12
MTetrahedron.h
fullMatrix::size1
int size1() const
Definition: fullMatrix.h:274
MTetrahedron::gammaShapeMeasure
virtual double gammaShapeMeasure()
Definition: MTetrahedron.cpp:94
z
const double z
Definition: GaussQuadratureQuad.cpp:56
picojson::copy
void copy(const std::string &s, Iter oi)
Definition: picojson.h:510
_myGetFaceRep
static void _myGetFaceRep(MTetrahedron *tet, int num, double *x, double *y, double *z, SVector3 *n, int numSubEdges)
Definition: MTetrahedron.cpp:214
qmTetrahedron::QMTET_GAMMA
@ QMTET_GAMMA
Definition: qualityMeasures.h:68
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
MTet4::circumcenter
void circumcenter(double *res)
Definition: meshGRegionDelaunayInsertion.h:75
MTetrahedronN::reverse
virtual void reverse()
Definition: MTetrahedron.cpp:366
MTetrahedronN::getNumEdgesRep
virtual int getNumEdgesRep(bool curved)
Definition: MTetrahedron.cpp:146
MVertex::y
double y() const
Definition: MVertex.h:61
pointsGenerators.h
MTetrahedron::etaShapeMeasure
virtual double etaShapeMeasure()
Definition: MTetrahedron.cpp:104
gmshGenerateMonomialsTetrahedron
fullMatrix< double > gmshGenerateMonomialsTetrahedron(int order, bool serendip)
Definition: pointsGenerators.cpp:317
MTetrahedron10::getEdgeRep
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MTetrahedron.cpp:186
SQU
#define SQU(a)
Definition: MTetrahedron.cpp:19
MVertex::distance
double distance(MVertex *const v)
Definition: MVertex.h:105
MTetrahedron::getFaceRep
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MTetrahedron.h:96
MTetrahedron::getFace
virtual MFace getFace(int num) const
Definition: MTetrahedron.h:88
MTetrahedronN::getEdgeRep
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MTetrahedron.cpp:195
MTetrahedron10::getNumFacesRep
virtual int getNumFacesRep(bool curved)
Definition: MTetrahedron.cpp:209
MVertex::x
double x() const
Definition: MVertex.h:60
SVector3::normalize
double normalize()
Definition: SVector3.h:38
IndicesReversed
std::vector< int > IndicesReversed
Definition: MHexahedron.h:605
BasisFactory.h