1 // Gmsh - Copyright (C) 1997-2021 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 <stdlib.h>
7 #include <cmath>
8 #include <limits>
9 #include "GmshConfig.h"
10 #include "GmshMessage.h"
11 #include "GModel.h"
12 #include "MElement.h"
13 #include "MPoint.h"
14 #include "MLine.h"
15 #include "MTriangle.h"
16 #include "MQuadrangle.h"
17 #include "MTetrahedron.h"
18 #include "MHexahedron.h"
19 #include "MPrism.h"
20 #include "MPyramid.h"
21 #include "MTrihedron.h"
22 #include "MElementCut.h"
23 #include "MSubElement.h"
24 #include "GEntity.h"
25 #include "StringUtils.h"
26 #include "Numeric.h"
27 #include "nodalBasis.h"
28 #include "CondNumBasis.h"
29 #include "Context.h"
30 #include "FuncSpaceData.h"
31 #include "bezierBasis.h"
32 #include "polynomialBasis.h"
33
34 #if defined(HAVE_MESH)
35 #include "qualityMeasuresJacobian.h"
36 #endif
37
38 #define SQU(a) ((a) * (a))
39
MElement(std::size_t num,int part)40 MElement::MElement(std::size_t num, int part) : _visible(1)
41 {
42 // we should make GModel a mandatory argument to the constructor
43 GModel *m = GModel::current();
44 if(num) {
45 _num = num;
46 m->setMaxElementNumber(_num);
47 }
48 else {
49 _num = m->incrementAndGetMaxElementNumber();
50 }
51 _partition = (short)part;
52 }
53
forceNum(std::size_t num)54 void MElement::forceNum(std::size_t num)
55 {
56 GModel *m = GModel::current();
57 _num = num;
58 m->setMaxElementNumber(_num);
59 }
60
getTolerance() const61 double MElement::getTolerance() const
62 {
63 return CTX::instance()->mesh.toleranceReferenceElement;
64 }
65
_getFaceInfo(const MFace & face,const MFace & other,int & sign,int & rot)66 bool MElement::_getFaceInfo(const MFace &face, const MFace &other, int &sign,
67 int &rot)
68 {
69 // Looks how is 'other' compared to 'face'. We suppose that 'face'
70 // is the reference.
71 // Return false if they are different.
72 // In case sign = 1, then rot = 1 if 'other' is equal to 'face' rotated
73 // one time according to the right hand side.
74 // In case sign = -1, then rot = 0 if 'other' is equal to reversed 'face'.
75 // In case sign = -1, then rot = 1 if 'other' is equal to reversed 'face'
76 // rotated one time according to the right hand side.
77 int N = face.getNumVertices();
78
79 sign = 0;
80 rot = -1;
81
82 if(face.getNumVertices() != other.getNumVertices()) return false;
83
84 sign = 1;
85 for(rot = 0; rot < N; ++rot) {
86 int i;
87 for(i = 0; i < N; ++i) {
88 if(other.getVertex(i) != face.getVertex((i + rot) % N)) break;
89 }
90 if(i == N) return true;
91 }
92
93 sign = -1;
94 for(rot = 0; rot < N; ++rot) {
95 int i;
96 for(i = 0; i < N; ++i) {
97 if(other.getVertex(i) != face.getVertex((N + rot - i) % N)) break;
98 }
99 if(i == N) return true;
100 }
101
102 sign = 0;
103 rot = -1;
104 return false;
105 }
106
_getEdgeRep(MVertex * v0,MVertex * v1,double * x,double * y,double * z,SVector3 * n,int faceIndex)107 void MElement::_getEdgeRep(MVertex *v0, MVertex *v1, double *x, double *y,
108 double *z, SVector3 *n, int faceIndex)
109 {
110 x[0] = v0->x();
111 y[0] = v0->y();
112 z[0] = v0->z();
113 x[1] = v1->x();
114 y[1] = v1->y();
115 z[1] = v1->z();
116 if(faceIndex >= 0) { n[0] = n[1] = getFace(faceIndex).normal(); }
117 else {
118 MEdge e(v0, v1);
119 n[0] = n[1] = e.normal();
120 }
121 }
122
123 #if defined(HAVE_VISUDEV)
_getFaceRepQuad(MVertex * v0,MVertex * v1,MVertex * v2,MVertex * v3,double * x,double * y,double * z,SVector3 * n)124 void MElement::_getFaceRepQuad(MVertex *v0, MVertex *v1, MVertex *v2,
125 MVertex *v3, double *x, double *y, double *z,
126 SVector3 *n)
127 {
128 x[0] = v0->x();
129 x[1] = v1->x();
130 x[2] = (x[0] + x[1] + v2->x() + v3->x()) / 4;
131 y[0] = v0->y();
132 y[1] = v1->y();
133 y[2] = (y[0] + y[1] + v2->y() + v3->y()) / 4;
134 z[0] = v0->z();
135 z[1] = v1->z();
136 z[2] = (z[0] + z[1] + v2->z() + v3->z()) / 4;
137 ;
138 SVector3 t1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
139 SVector3 t2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
140 SVector3 normal = crossprod(t1, t2);
141 normal.normalize();
142 for(int i = 0; i < 3; i++) n[i] = normal;
143 }
144 #endif
145
_getFaceRep(MVertex * v0,MVertex * v1,MVertex * v2,double * x,double * y,double * z,SVector3 * n)146 void MElement::_getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2, double *x,
147 double *y, double *z, SVector3 *n)
148 {
149 x[0] = v0->x();
150 x[1] = v1->x();
151 x[2] = v2->x();
152 y[0] = v0->y();
153 y[1] = v1->y();
154 y[2] = v2->y();
155 z[0] = v0->z();
156 z[1] = v1->z();
157 z[2] = v2->z();
158 SVector3 t1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
159 SVector3 t2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
160 SVector3 normal = crossprod(t1, t2);
161 normal.normalize();
162 for(int i = 0; i < 3; i++) n[i] = normal;
163 }
164
getVisibility() const165 char MElement::getVisibility() const
166 {
167 if(CTX::instance()->hideUnselected && _visible < 2) return false;
168 return _visible;
169 }
170
getHighOrderEdge(int num,int sign)171 MEdgeN MElement::getHighOrderEdge(int num, int sign)
172 {
173 const int order = getPolynomialOrder();
174 std::vector<MVertex *> vertices(static_cast<std::size_t>(order) + 1);
175 vertices[0] = getVertex(numEdge2numVertex(num, sign > 0 ? 0 : 1));
176 vertices[1] = getVertex(numEdge2numVertex(num, sign > 0 ? 1 : 0));
177 const int start = (int)getNumPrimaryVertices() + num * (order - 1);
178 const int end = (int)getNumPrimaryVertices() + (num + 1) * (order - 1);
179 int k = 1;
180 if(sign > 0) {
181 for(int i = start; i < end; ++i) { vertices[++k] = getVertex(i); }
182 }
183 else {
184 for(int i = end - 1; i >= start; --i) { vertices[++k] = getVertex(i); }
185 }
186 return MEdgeN(vertices);
187 }
188
getEdgeInfo(const MEdge & edge,int & ithEdge,int & sign) const189 bool MElement::getEdgeInfo(const MEdge &edge, int &ithEdge, int &sign) const
190 {
191 for(ithEdge = 0; ithEdge < getNumEdges(); ithEdge++) {
192 const MVertex *v0 = getVertex(numEdge2numVertex(ithEdge, 0));
193 const MVertex *v1 = getVertex(numEdge2numVertex(ithEdge, 1));
194 if(v0 == edge.getVertex(0) && v1 == edge.getVertex(1)) {
195 sign = 1;
196 return true;
197 }
198 if(v1 == edge.getVertex(0) && v0 == edge.getVertex(1)) {
199 sign = -1;
200 return true;
201 }
202 }
203 Msg::Error("Could not get edge information for element %lu", getNum());
204 return false;
205 }
206
getHighOrderFace(int num,int sign,int rot)207 MFaceN MElement::getHighOrderFace(int num, int sign, int rot)
208 {
209 if(getDim() < 2 || getDim() > 3) {
210 Msg::Error("Wrong dimension for getHighOrderFace");
211 return MFaceN();
212 }
213
214 if(getDim() == 2) {
215 std::vector<MVertex *> vertices(getNumVertices());
216 getVertices(vertices);
217 return MFaceN(getType(), getPolynomialOrder(), vertices);
218 }
219
220 const nodalBasis *fs = getFunctionSpace();
221 int id = fs->getClosureId(num, sign, rot);
222 const std::vector<int> &closure = fs->getClosure(id);
223
224 std::vector<MVertex *> vertices(closure.size());
225 for(std::size_t i = 0; i < closure.size(); ++i) {
226 vertices[i] = getVertex(closure[i]);
227 }
228
229 static int type2numTriFaces[9] = {0, 0, 0, 1, 0, 4, 4, 2, 0};
230 int typeFace = TYPE_TRI;
231 if(num >= type2numTriFaces[getType()]) typeFace = TYPE_QUA;
232
233 return MFaceN(typeFace, getPolynomialOrder(), vertices);
234 }
235
minEdge()236 double MElement::minEdge()
237 {
238 double m = 1.e25;
239 for(int i = 0; i < getNumEdges(); i++) {
240 MEdge e = getEdge(i);
241 m = std::min(m, e.getVertex(0)->distance(e.getVertex(1)));
242 }
243 return m;
244 }
245
maxEdge()246 double MElement::maxEdge()
247 {
248 double m = 0.;
249 for(int i = 0; i < getNumEdges(); i++) {
250 MEdge e = getEdge(i);
251 m = std::max(m, e.getVertex(0)->distance(e.getVertex(1)));
252 }
253 return m;
254 }
255
maxDistToStraight() const256 double MElement::maxDistToStraight() const
257 {
258 const nodalBasis *lagBasis = getFunctionSpace();
259 const fullMatrix<double> &uvw = lagBasis->points;
260 const int &nV = uvw.size1();
261 const int &dim = uvw.size2();
262 const nodalBasis *lagBasis1 = getFunctionSpace(1);
263 const int &nV1 = lagBasis1->points.size1();
264 std::vector<SPoint3> xyz1(nV1);
265 for(int iV = 0; iV < nV1; ++iV) xyz1[iV] = getVertex(iV)->point();
266 double maxdx = 0.;
267 for(int iV = nV1; iV < nV; ++iV) {
268 double f[256];
269 lagBasis1->f(uvw(iV, 0), (dim > 1) ? uvw(iV, 1) : 0.,
270 (dim > 2) ? uvw(iV, 2) : 0., f);
271 SPoint3 xyzS(0., 0., 0.);
272 for(int iSF = 0; iSF < nV1; ++iSF) xyzS += xyz1[iSF] * f[iSF];
273 SVector3 vec(xyzS, getVertex(iV)->point());
274 double dx = vec.norm();
275 if(dx > maxdx) maxdx = dx;
276 }
277 return maxdx;
278 }
279
minIsotropyMeasure(bool knownValid,bool reversedOK)280 double MElement::minIsotropyMeasure(bool knownValid, bool reversedOK)
281 {
282 #if defined(HAVE_MESH)
283 return jacobianBasedQuality::minICNMeasure(this, knownValid, reversedOK);
284 #else
285 return 0.;
286 #endif
287 }
288
minScaledJacobian(bool knownValid,bool reversedOK)289 double MElement::minScaledJacobian(bool knownValid, bool reversedOK)
290 {
291 #if defined(HAVE_MESH)
292 return jacobianBasedQuality::minIGEMeasure(this, knownValid, reversedOK);
293 #else
294 return 0.;
295 #endif
296 }
297
scaledJacRange(double & jmin,double & jmax,GEntity * ge) const298 void MElement::scaledJacRange(double &jmin, double &jmax, GEntity *ge) const
299 {
300 jmin = jmax = 1.0;
301 #if defined(HAVE_MESH)
302 const JacobianBasis *jac = getJacobianFuncSpace();
303 fullMatrix<double> nodesXYZ(jac->getNumMapNodes(), 3);
304 getNodesCoord(nodesXYZ);
305 fullVector<double> SJi(jac->getNumSamplingPnts());
306 jac->getScaledJacobian(nodesXYZ, SJi);
307 if(ge && (ge->dim() == 2) && ge->haveParametrization()) {
308 // If parametrized surface entity provided...
309 SVector3 geoNorm(0., 0., 0.);
310 // ... correct Jacobian sign with geometrical normal
311 for(int i = 0; i < jac->getNumPrimMapNodes(); i++) {
312 const MVertex *vert = getVertex(i);
313 if(vert->onWhat() == ge) {
314 double u, v;
315 vert->getParameter(0, u);
316 vert->getParameter(1, v);
317 geoNorm += ((GFace *)ge)->normal(SPoint2(u, v));
318 }
319 }
320 if(geoNorm.normSq() == 0.) {
321 // If no vertex on surface or average is zero, take normal at barycenter
322 SPoint2 param = ((GFace *)ge)->parFromPoint(barycenter(true), false);
323 geoNorm = ((GFace *)ge)->normal(param);
324 }
325 fullMatrix<double> elNorm(1, 3);
326 jac->getPrimNormal2D(nodesXYZ, elNorm);
327 const double scal = geoNorm(0) * elNorm(0, 0) + geoNorm(1) * elNorm(0, 1) +
328 geoNorm(2) * elNorm(0, 2);
329 if(scal < 0.) SJi.scale(-1.);
330 }
331 bezierCoeff Bi(jac->getFuncSpaceData(), SJi);
332 jmin = *std::min_element(Bi.getDataPtr(), Bi.getDataPtr() + Bi.getNumCoeff());
333 jmax = *std::max_element(Bi.getDataPtr(), Bi.getDataPtr() + Bi.getNumCoeff());
334 #endif
335 }
336
idealJacRange(double & jmin,double & jmax,GEntity * ge)337 void MElement::idealJacRange(double &jmin, double &jmax, GEntity *ge)
338 {
339 jmin = jmax = 1.0;
340 #if defined(HAVE_MESH)
341 const JacobianBasis *jac = getJacobianFuncSpace();
342 fullMatrix<double> nodesXYZ(jac->getNumMapNodes(), 3);
343 getNodesCoord(nodesXYZ);
344 fullVector<double> iJi(jac->getNumSamplingPnts());
345 jac->getSignedIdealJacobian(nodesXYZ, iJi);
346 const int nEd = getNumEdges(), dim = getDim();
347 double sumEdLength = 0.;
348 for(int iEd = 0; iEd < nEd; iEd++) sumEdLength += getEdge(iEd).length();
349 const double invMeanEdLength = double(nEd) / sumEdLength;
350 if(sumEdLength == 0.) {
351 jmin = 0.;
352 jmax = 0.;
353 return;
354 }
355 double scale =
356 (dim == 1.) ? invMeanEdLength :
357 (dim == 2.) ? invMeanEdLength * invMeanEdLength :
358 invMeanEdLength * invMeanEdLength * invMeanEdLength;
359 if(ge && (ge->dim() == 2) && ge->haveParametrization()) {
360 // If parametrized surface entity provided...
361 SVector3 geoNorm(0., 0., 0.);
362 // ... correct Jacobian sign with geometrical normal
363 for(int i = 0; i < jac->getNumPrimMapNodes(); i++) {
364 const MVertex *vert = getVertex(i);
365 if(vert->onWhat() == ge) {
366 double u, v;
367 vert->getParameter(0, u);
368 vert->getParameter(1, v);
369 geoNorm += ((GFace *)ge)->normal(SPoint2(u, v));
370 }
371 }
372 if(geoNorm.normSq() == 0.) {
373 // If no vertex on surface or average is zero, take normal at barycenter
374 SPoint2 param = ((GFace *)ge)->parFromPoint(barycenter(true), false);
375 geoNorm = ((GFace *)ge)->normal(param);
376 }
377 fullMatrix<double> elNorm(1, 3);
378 jac->getPrimNormal2D(nodesXYZ, elNorm, true);
379 const double dp = geoNorm(0) * elNorm(0, 0) + geoNorm(1) * elNorm(0, 1) +
380 geoNorm(2) * elNorm(0, 2);
381 if(dp < 0.) scale = -scale;
382 }
383 bezierCoeff Bi(jac->getFuncSpaceData(), iJi);
384 jmin = *std::min_element(Bi.getDataPtr(), Bi.getDataPtr() + Bi.getNumCoeff());
385 jmax = *std::max_element(Bi.getDataPtr(), Bi.getDataPtr() + Bi.getNumCoeff());
386 #endif
387 }
388
signedInvCondNumRange(double & iCNMin,double & iCNMax,GEntity * ge)389 void MElement::signedInvCondNumRange(double &iCNMin, double &iCNMax,
390 GEntity *ge)
391 {
392 iCNMin = iCNMax = 1.0;
393 #if defined(HAVE_MESH)
394 const CondNumBasis *cnb = BasisFactory::getCondNumBasis(getTypeForMSH());
395 const int numCNNodes = cnb->getNumCondNumNodes();
396 fullMatrix<double> nodesXYZ(cnb->getNumMapNodes(), 3), normals;
397 getNodesCoord(nodesXYZ);
398 if(getDim() == 2.) {
399 SVector3 nVec = getFace(0).normal();
400 normals.resize(1, 3);
401 normals(0, 0) = nVec[0];
402 normals(0, 1) = nVec[1];
403 normals(0, 2) = nVec[2];
404 }
405 if(ge && (ge->dim() == 2) && ge->haveParametrization()) {
406 // If parametrized surface entity provided...
407 SVector3 geoNorm(0., 0., 0.);
408 // ... correct Jacobian sign with geometrical normal
409 for(std::size_t i = 0; i < getNumPrimaryVertices(); i++) {
410 const MVertex *vert = getVertex(i);
411 if(vert->onWhat() == ge) {
412 double u, v;
413 vert->getParameter(0, u);
414 vert->getParameter(1, v);
415 geoNorm += ((GFace *)ge)->normal(SPoint2(u, v));
416 }
417 }
418 if(geoNorm.normSq() == 0.) {
419 // If no vertex on surface or average is zero, take normal at barycenter
420 SPoint2 param = ((GFace *)ge)->parFromPoint(barycenter(true), false);
421 geoNorm = ((GFace *)ge)->normal(param);
422 }
423 const double dp = geoNorm(0) * normals(0, 0) + geoNorm(1) * normals(0, 1) +
424 geoNorm(2) * normals(0, 2);
425 if(dp < 0.) {
426 normals(0, 0) = -normals(0, 0);
427 normals(0, 1) = -normals(0, 1);
428 normals(0, 2) = -normals(0, 2);
429 }
430 }
431 fullVector<double> invCondNum(numCNNodes);
432 cnb->getSignedInvCondNum(nodesXYZ, normals, invCondNum);
433 iCNMin = *std::min_element(invCondNum.getDataPtr(),
434 invCondNum.getDataPtr() + numCNNodes);
435 iCNMax = *std::max_element(invCondNum.getDataPtr(),
436 invCondNum.getDataPtr() + numCNNodes);
437 #endif
438 }
439
signedInvGradErrorRange(double & minSIGE,double & maxSIGE)440 void MElement::signedInvGradErrorRange(double &minSIGE, double &maxSIGE)
441 {
442 #if defined(HAVE_MESH)
443 jacobianBasedQuality::sampleIGEMeasure(this, getPolynomialOrder(), minSIGE,
444 maxSIGE);
445 #endif
446 }
447
getNode(int num,double & u,double & v,double & w) const448 void MElement::getNode(int num, double &u, double &v, double &w) const
449 {
450 // Should we always do this instead of using lookup table for linear elements?
451 const nodalBasis *nb = getFunctionSpace();
452 const fullMatrix<double> &refpnts = nb->getReferenceNodes();
453 u = refpnts(num, 0);
454 v = getDim() > 1 ? refpnts(num, 1) : 0;
455 w = getDim() > 2 ? refpnts(num, 2) : 0;
456 }
457
getShapeFunctions(double u,double v,double w,double s[],int o) const458 void MElement::getShapeFunctions(double u, double v, double w, double s[],
459 int o) const
460 {
461 const nodalBasis *fs = getFunctionSpace(o);
462 if(fs)
463 fs->f(u, v, w, s);
464 else
465 Msg::Error("Function space not implemented for this type of element");
466 }
467
getGradShapeFunctions(double u,double v,double w,double s[][3],int o) const468 void MElement::getGradShapeFunctions(double u, double v, double w,
469 double s[][3], int o) const
470 {
471 const nodalBasis *fs = getFunctionSpace(o);
472 if(fs)
473 fs->df(u, v, w, s);
474 else
475 Msg::Error("Function space not implemented for this type of element");
476 }
477
getHessShapeFunctions(double u,double v,double w,double s[][3][3],int o) const478 void MElement::getHessShapeFunctions(double u, double v, double w,
479 double s[][3][3], int o) const
480 {
481 const nodalBasis *fs = getFunctionSpace(o);
482 if(fs)
483 fs->ddf(u, v, w, s);
484 else
485 Msg::Error("Function space not implemented for this type of element");
486 }
487
getThirdDerivativeShapeFunctions(double u,double v,double w,double s[][3][3][3],int o) const488 void MElement::getThirdDerivativeShapeFunctions(double u, double v, double w,
489 double s[][3][3][3],
490 int o) const
491 {
492 const nodalBasis *fs = getFunctionSpace(o);
493 if(fs)
494 fs->dddf(u, v, w, s);
495 else
496 Msg::Error("Function space not implemented for this type of element");
497 }
498
barycenter_infty() const499 SPoint3 MElement::barycenter_infty() const
500 {
501 double xmin = getVertex(0)->x();
502 double xmax = xmin;
503 double ymin = getVertex(0)->y();
504 double ymax = ymin;
505 double zmin = getVertex(0)->z();
506 double zmax = zmin;
507 int n = getNumVertices();
508 for(int i = 0; i < n; i++) {
509 const MVertex *v = getVertex(i);
510 xmin = std::min(xmin, v->x());
511 xmax = std::max(xmax, v->x());
512 ymin = std::min(ymin, v->y());
513 ymax = std::max(ymax, v->y());
514 zmin = std::min(zmin, v->z());
515 zmax = std::max(zmax, v->z());
516 }
517 return SPoint3(0.5 * (xmin + xmax), 0.5 * (ymin + ymax), 0.5 * (zmin + zmax));
518 }
519
barycenter(bool primary) const520 SPoint3 MElement::barycenter(bool primary) const
521 {
522 SPoint3 p(0., 0., 0.);
523 std::size_t n = primary ? getNumPrimaryVertices() : getNumVertices();
524 for(std::size_t i = 0; i < n; i++) {
525 const MVertex *v = getVertex(i);
526 p[0] += v->x();
527 p[1] += v->y();
528 p[2] += v->z();
529 }
530 p[0] /= (double)n;
531 p[1] /= (double)n;
532 p[2] /= (double)n;
533 return p;
534 }
535
fastBarycenter(bool primary) const536 SPoint3 MElement::fastBarycenter(bool primary) const
537 {
538 SPoint3 p(0., 0., 0.);
539 std::size_t n = primary ? getNumPrimaryVertices() : getNumVertices();
540 for(std::size_t i = 0; i < n; i++) {
541 const MVertex *v = getVertex(i);
542 p[0] += v->x();
543 p[1] += v->y();
544 p[2] += v->z();
545 }
546
547 return p;
548 }
549
barycenterUVW() const550 SPoint3 MElement::barycenterUVW() const
551 {
552 SPoint3 p(0., 0., 0.);
553 int n = getNumVertices();
554 for(int i = 0; i < n; i++) {
555 double x, y, z;
556 getNode(i, x, y, z);
557 p[0] += x;
558 p[1] += y;
559 p[2] += z;
560 }
561 p[0] /= (double)n;
562 p[1] /= (double)n;
563 p[2] /= (double)n;
564 return p;
565 }
566
getVolume()567 double MElement::getVolume()
568 {
569 int npts;
570 IntPt *pts;
571 getIntegrationPoints(getDim() * (getPolynomialOrder() - 1), &npts, &pts);
572 double vol = 0.;
573 for(int i = 0; i < npts; i++) {
574 vol += getJacobianDeterminant(pts[i].pt[0], pts[i].pt[1], pts[i].pt[2]) *
575 pts[i].weight;
576 }
577 return vol;
578 }
579
getVolumeSign()580 int MElement::getVolumeSign()
581 {
582 double v = getVolume();
583 if(v < 0.)
584 return -1;
585 else if(v > 0.)
586 return 1;
587 else
588 return 0;
589 }
590
setVolumePositive()591 bool MElement::setVolumePositive()
592 {
593 if(getDim() < 3) return true;
594 int s = getVolumeSign();
595 if(s < 0) reverse();
596 if(!s) return false;
597 return true;
598 }
599
getValidity()600 int MElement::getValidity()
601 {
602 #if defined(HAVE_MESH)
603 double jmin, jmax;
604 jacobianBasedQuality::minMaxJacobianDeterminant(this, jmin, jmax);
605 if(jmin > .0) return 1; // valid
606 if(jmax >= .0) return 0; // invalid
607 // Here, jmax < 0 (and jmin < 0). The element validity is quite indeterminate.
608 // It can be valid but with a wrong numbering of the nodes,
609 // or it can be invalid, i.e. with nodes that are incorrectly located.
610 return -1;
611 #else
612 return 0;
613 #endif
614 }
615
getInfoString(bool multline)616 std::string MElement::getInfoString(bool multline)
617 {
618 std::ostringstream sstream;
619 sstream.precision(12);
620
621 sstream << "Element " << getNum() << ":";
622 if(multline) sstream << "\n";
623
624 const char *name;
625 MElement::getInfoMSH(getTypeForMSH(), &name);
626 sstream << " " << name << " (MSH type " << getTypeForMSH() << ", dimension "
627 << getDim() << ", order " << getPolynomialOrder() << ", partition "
628 << getPartition() << ")";
629 if(multline) sstream << "\n";
630
631 sstream << " Nodes:";
632 for(std::size_t i = 0; i < getNumVertices(); i++)
633 sstream << " " << getVertex(i)->getNum();
634 if(multline) sstream << "\n";
635
636 sstream << " Volume: " << getVolume() << "\n";
637
638 SPoint3 pt = barycenter();
639 sstream << " Barycenter: (" << pt[0] << ", " << pt[1] << ", " << pt[2] << ")";
640 if(multline) sstream << "\n";
641
642 sstream << " Edge length: "
643 << "min = " << minEdge() << " "
644 << "max = " << maxEdge();
645 if(multline) sstream << "\n";
646
647 sstream << " Quality: "
648 << "gamma = " << gammaShapeMeasure();
649 if(multline) sstream << "\n";
650
651 double sICNMin, sICNMax;
652 signedInvCondNumRange(sICNMin, sICNMax);
653 sstream << " SICN range: " << sICNMin << " " << sICNMax;
654 if(multline) sstream << "\n";
655
656 double sIGEMin, sIGEMax;
657 signedInvGradErrorRange(sIGEMin, sIGEMax);
658 sstream << " SIGE range: " << sIGEMin << " " << sIGEMax;
659 if(multline) sstream << "\n";
660
661 sstream << " Inner / outer radius: " << getInnerRadius() << " / "
662 << getOuterRadius();
663 return sstream.str();
664 }
665
getFunctionSpace(int order,bool serendip) const666 const nodalBasis *MElement::getFunctionSpace(int order, bool serendip) const
667 {
668 if(order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
669 int type = ElementType::getType(getType(), order, serendip);
670 return type ? BasisFactory::getNodalBasis(type) : nullptr;
671 }
672
getFuncSpaceData(int order,bool serendip) const673 const FuncSpaceData MElement::getFuncSpaceData(int order, bool serendip) const
674 {
675 if(order == -1) return FuncSpaceData(this);
676 return FuncSpaceData(this, order, serendip);
677 }
678
getJacobianFuncSpace(int orderElement) const679 const JacobianBasis *MElement::getJacobianFuncSpace(int orderElement) const
680 {
681 if(orderElement == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
682 int tag = ElementType::getType(getType(), orderElement);
683 return tag ? BasisFactory::getJacobianBasis(tag) : nullptr;
684 }
685
getJacobianFuncSpaceData(int orderElement) const686 const FuncSpaceData MElement::getJacobianFuncSpaceData(int orderElement) const
687 {
688 if(orderElement == -1) orderElement = getPolynomialOrder();
689 int orderJac = JacobianBasis::jacobianOrder(this->getType(), orderElement);
690 return FuncSpaceData(this, orderJac, false);
691 }
692
_computeDeterminantAndRegularize(const MElement * ele,double jac[3][3])693 static double _computeDeterminantAndRegularize(const MElement *ele,
694 double jac[3][3])
695 {
696 double dJ = 0;
697
698 switch(ele->getDim()) {
699 case 0: {
700 dJ = 1.0;
701 jac[0][0] = jac[1][1] = jac[2][2] = 1.0;
702 jac[0][1] = jac[1][0] = jac[2][0] = 0.0;
703 jac[0][2] = jac[1][2] = jac[2][1] = 0.0;
704 break;
705 }
706 case 1: {
707 dJ = sqrt(SQU(jac[0][0]) + SQU(jac[0][1]) + SQU(jac[0][2]));
708
709 // regularize matrix
710 double a[3], b[3], c[3];
711 a[0] = jac[0][0];
712 a[1] = jac[0][1];
713 a[2] = jac[0][2];
714 if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
715 (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))) {
716 b[0] = a[1];
717 b[1] = -a[0];
718 b[2] = 0.;
719 }
720 else {
721 b[0] = 0.;
722 b[1] = a[2];
723 b[2] = -a[1];
724 }
725 norme(b);
726 prodve(a, b, c);
727 norme(c);
728 jac[1][0] = b[0];
729 jac[1][1] = b[1];
730 jac[1][2] = b[2];
731 jac[2][0] = c[0];
732 jac[2][1] = c[1];
733 jac[2][2] = c[2];
734 break;
735 }
736 case 2: {
737 dJ = sqrt(SQU(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
738 SQU(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
739 SQU(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
740
741 // regularize matrix
742 double a[3], b[3], c[3];
743 a[0] = jac[0][0];
744 a[1] = jac[0][1];
745 a[2] = jac[0][2];
746 b[0] = jac[1][0];
747 b[1] = jac[1][1];
748 b[2] = jac[1][2];
749 prodve(a, b, c);
750 norme(c);
751 jac[2][0] = c[0];
752 jac[2][1] = c[1];
753 jac[2][2] = c[2];
754 break;
755 }
756 case 3: {
757 dJ =
758 (jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
759 jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
760 jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2]);
761 break;
762 }
763 }
764 return dJ;
765 }
766
_computeDeterminantAndRegularize(const MElement * ele,double * jac)767 static double _computeDeterminantAndRegularize(const MElement *ele, double *jac)
768 {
769 double dJ = 0;
770
771 // 'jac' is a row-major order array :
772 //
773 // |0 1 2|
774 // |3 4 5|
775 // |6 7 8|
776
777 switch(ele->getDim()) {
778 case 0: {
779 dJ = 1.0;
780 jac[0] = jac[4] = jac[8] = 1.0;
781 jac[1] = jac[2] = jac[3] = 0.0;
782 jac[5] = jac[6] = jac[7] = 0.0;
783 break;
784 }
785 case 1: {
786 dJ = sqrt(SQU(jac[0]) + SQU(jac[1]) + SQU(jac[2]));
787
788 // regularize matrix
789 double a[3], b[3], c[3];
790 a[0] = jac[0];
791 a[1] = jac[1];
792 a[2] = jac[2];
793 if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
794 (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))) {
795 b[0] = a[1];
796 b[1] = -a[0];
797 b[2] = 0.;
798 }
799 else {
800 b[0] = 0.;
801 b[1] = a[2];
802 b[2] = -a[1];
803 }
804 norme(b);
805 prodve(a, b, c);
806 norme(c);
807 jac[3] = b[0];
808 jac[4] = b[1];
809 jac[5] = b[2];
810 jac[6] = c[0];
811 jac[7] = c[1];
812 jac[8] = c[2];
813 break;
814 }
815 case 2: {
816 dJ = sqrt(SQU(jac[0] * jac[4] - jac[1] * jac[3]) +
817 SQU(jac[2] * jac[3] - jac[0] * jac[5]) +
818 SQU(jac[1] * jac[5] - jac[2] * jac[4]));
819
820 // regularize matrix
821 double a[3], b[3], c[3];
822 a[0] = jac[0];
823 a[1] = jac[1];
824 a[2] = jac[2];
825 b[0] = jac[3];
826 b[1] = jac[4];
827 b[2] = jac[5];
828 prodve(a, b, c);
829 norme(c);
830 jac[6] = c[0];
831 jac[7] = c[1];
832 jac[8] = c[2];
833 break;
834 }
835 case 3: {
836 dJ = (jac[0] * jac[4] * jac[8] + jac[2] * jac[3] * jac[7] +
837 jac[1] * jac[5] * jac[6] - jac[2] * jac[4] * jac[6] -
838 jac[0] * jac[5] * jac[7] - jac[1] * jac[3] * jac[8]);
839 break;
840 }
841 }
842 return dJ;
843 }
844
getJacobian(double u,double v,double w,double jac[3][3]) const845 double MElement::getJacobian(double u, double v, double w,
846 double jac[3][3]) const
847 {
848 jac[0][0] = jac[0][1] = jac[0][2] = 0.;
849 jac[1][0] = jac[1][1] = jac[1][2] = 0.;
850 jac[2][0] = jac[2][1] = jac[2][2] = 0.;
851 if(getDim() > 3) return 0;
852
853 double gsf[1256][3];
854 getGradShapeFunctions(u, v, w, gsf);
855 for(std::size_t i = 0; i < getNumShapeFunctions(); i++) {
856 const MVertex *ver = getShapeFunctionNode(i);
857 double *gg = gsf[i];
858 for(int j = 0; j < getDim(); j++) {
859 jac[j][0] += ver->x() * gg[j];
860 jac[j][1] += ver->y() * gg[j];
861 jac[j][2] += ver->z() * gg[j];
862 }
863 }
864
865 return _computeDeterminantAndRegularize(this, jac);
866 }
867
getJacobian(const fullMatrix<double> & gsf,double jac[3][3]) const868 double MElement::getJacobian(const fullMatrix<double> &gsf,
869 double jac[3][3]) const
870 {
871 jac[0][0] = jac[0][1] = jac[0][2] = 0.;
872 jac[1][0] = jac[1][1] = jac[1][2] = 0.;
873 jac[2][0] = jac[2][1] = jac[2][2] = 0.;
874 if(gsf.size2() > 3) return 0;
875
876 const std::size_t numShapeFunctions = getNumShapeFunctions();
877 for(std::size_t i = 0; i < numShapeFunctions; i++) {
878 const MVertex *v = getShapeFunctionNode(i);
879 for(int j = 0; j < gsf.size2(); j++) {
880 jac[j][0] += v->x() * gsf(i, j);
881 jac[j][1] += v->y() * gsf(i, j);
882 jac[j][2] += v->z() * gsf(i, j);
883 }
884 }
885 return _computeDeterminantAndRegularize(this, jac);
886 }
887
getJacobian(const std::vector<SVector3> & gsf,double jac[3][3]) const888 double MElement::getJacobian(const std::vector<SVector3> &gsf,
889 double jac[3][3]) const
890 {
891 jac[0][0] = jac[0][1] = jac[0][2] = 0.;
892 jac[1][0] = jac[1][1] = jac[1][2] = 0.;
893 jac[2][0] = jac[2][1] = jac[2][2] = 0.;
894
895 const int numShapeFunctions = getNumVertices();
896 for(int i = 0; i < numShapeFunctions; i++) {
897 const MVertex *v = getShapeFunctionNode(i);
898 for(int j = 0; j < 3; j++) {
899 const double mult = gsf[i][j];
900 jac[j][0] += v->x() * mult;
901 jac[j][1] += v->y() * mult;
902 jac[j][2] += v->z() * mult;
903 }
904 }
905 return _computeDeterminantAndRegularize(this, jac);
906 }
907
getJacobian(const std::vector<SVector3> & gsf,double * jac) const908 double MElement::getJacobian(const std::vector<SVector3> &gsf,
909 double *jac) const
910 {
911 for(std::size_t i = 0; i < 9; i++) { jac[i] = 0.; }
912
913 const int numShapeFunctions = getNumVertices();
914 for(int i = 0; i < numShapeFunctions; i++) {
915 const MVertex *v = getShapeFunctionNode(i);
916 for(int j = 0; j < 3; j++) {
917 const double mult = gsf[i][j];
918 jac[3 * j + 0] += v->x() * mult;
919 jac[3 * j + 1] += v->y() * mult;
920 jac[3 * j + 2] += v->z() * mult;
921 }
922 }
923 return _computeDeterminantAndRegularize(this, jac);
924 }
925
getJacobian(double u,double v,double w,fullMatrix<double> & j) const926 double MElement::getJacobian(double u, double v, double w,
927 fullMatrix<double> &j) const
928 {
929 double JAC[3][3];
930 const double detJ = getJacobian(u, v, w, JAC);
931 for(int i = 0; i < 3; i++) {
932 j(i, 0) = JAC[i][0];
933 j(i, 1) = JAC[i][1];
934 j(i, 2) = JAC[i][2];
935 }
936 return detJ;
937 }
938
getPrimaryJacobian(double u,double v,double w,double jac[3][3]) const939 double MElement::getPrimaryJacobian(double u, double v, double w,
940 double jac[3][3]) const
941 {
942 jac[0][0] = jac[0][1] = jac[0][2] = 0.;
943 jac[1][0] = jac[1][1] = jac[1][2] = 0.;
944 jac[2][0] = jac[2][1] = jac[2][2] = 0.;
945
946 double gsf[1256][3];
947 getGradShapeFunctions(u, v, w, gsf, 1);
948 for(std::size_t i = 0; i < getNumPrimaryShapeFunctions(); i++) {
949 const MVertex *v = getShapeFunctionNode(i);
950 double *gg = gsf[i];
951 for(std::size_t j = 0; j < 3; j++) {
952 jac[j][0] += v->x() * gg[j];
953 jac[j][1] += v->y() * gg[j];
954 jac[j][2] += v->z() * gg[j];
955 }
956 }
957
958 return _computeDeterminantAndRegularize(this, jac);
959 }
960
getSignedJacobian(fullVector<double> & jacobian,int o) const961 void MElement::getSignedJacobian(fullVector<double> &jacobian, int o) const
962 {
963 const int numNodes = getNumVertices();
964 fullMatrix<double> nodesXYZ(numNodes, 3);
965 getNodesCoord(nodesXYZ);
966 getJacobianFuncSpace(o)->getSignedJacobian(nodesXYZ, jacobian);
967 }
968
getNodesCoord(fullMatrix<double> & nodesXYZ) const969 void MElement::getNodesCoord(fullMatrix<double> &nodesXYZ) const
970 {
971 const int numNodes = getNumVertices();
972 for(int i = 0; i < numNodes; i++) {
973 const MVertex *v = getShapeFunctionNode(i);
974 nodesXYZ(i, 0) = v->x();
975 nodesXYZ(i, 1) = v->y();
976 nodesXYZ(i, 2) = v->z();
977 }
978 }
979
getNodesCoordNonSerendip(fullMatrix<double> & nodesXYZ) const980 void MElement::getNodesCoordNonSerendip(fullMatrix<double> &nodesXYZ) const
981 {
982 int tag = ElementType::getType(getType(), getPolynomialOrder(), false);
983 const nodalBasis *basis = BasisFactory::getNodalBasis(tag);
984 const fullMatrix<double> nodesUVW = basis->getReferenceNodes();
985 nodesXYZ.resize(nodesUVW.size1(), 3, false);
986
987 getNodesCoord(nodesXYZ);
988
989 double xyz[3];
990 for(int i = getNumVertices(); i < nodesUVW.size1(); ++i) {
991 pnt(nodesUVW(i, 0), nodesUVW(i, 1), nodesUVW(i, 2), xyz);
992 nodesXYZ(i, 0) = xyz[0];
993 nodesXYZ(i, 1) = xyz[1];
994 nodesXYZ(i, 2) = xyz[2];
995 }
996 }
997
getBezierVerticesCoord() const998 bezierCoeff *MElement::getBezierVerticesCoord() const
999 {
1000 const bezierBasis *basis;
1001 basis = BasisFactory::getBezierBasis(getType(), getPolynomialOrder());
1002 const fullMatrix<double> pntUVW =
1003 basis->getSamplingPointsToComputeBezierCoeff();
1004
1005 fullMatrix<double> pntXYZ(pntUVW.size1(), 3);
1006 double xyz[3];
1007 double uvw[3];
1008 for(int i = 0; i < pntUVW.size1(); ++i) {
1009 uvw[0] = uvw[1] = uvw[2] = 0.;
1010 for(int j = 0; j < pntUVW.size2(); ++j) {
1011 uvw[j] = pntUVW(i, j);
1012 }
1013 pnt(uvw[0], uvw[1], uvw[2], xyz);
1014 pntXYZ(i, 0) = xyz[0];
1015 pntXYZ(i, 1) = xyz[1];
1016 pntXYZ(i, 2) = xyz[2];
1017 }
1018
1019 return new bezierCoeff(getFuncSpaceData(getPolynomialOrder(), false), pntXYZ);
1020 }
1021
getEigenvaluesMetric(double u,double v,double w,double values[3]) const1022 double MElement::getEigenvaluesMetric(double u, double v, double w,
1023 double values[3]) const
1024 {
1025 double jac[3][3];
1026 getJacobian(u, v, w, jac);
1027 GradientBasis::mapFromIdealElement(getType(), jac);
1028
1029 switch(getDim()) {
1030 case 1:
1031 values[0] = 0;
1032 values[1] = -1;
1033 values[2] = -1;
1034 for(int d = 0; d < 3; ++d) values[0] += jac[d][0] * jac[d][0];
1035 return 1;
1036
1037 case 2: {
1038 fullMatrix<double> metric(2, 2);
1039 for(int i = 0; i < 2; ++i) {
1040 for(int j = 0; j < 2; ++j) {
1041 for(int d = 0; d < 3; ++d) metric(i, j) += jac[d][i] * jac[d][j];
1042 }
1043 }
1044 fullVector<double> valReal(values, 2), valImag(2);
1045 fullMatrix<double> vecLeft(2, 2), vecRight(2, 2);
1046 metric.eig(valReal, valImag, vecLeft, vecRight, true);
1047 values[2] = -1;
1048 return std::sqrt(valReal(0) / valReal(1));
1049 }
1050
1051 case 3: {
1052 fullMatrix<double> metric(3, 3);
1053 for(int i = 0; i < 3; ++i) {
1054 for(int j = 0; j < 3; ++j) {
1055 for(int d = 0; d < 3; ++d) metric(i, j) += jac[d][i] * jac[d][j];
1056 }
1057 }
1058
1059 fullVector<double> valReal(values, 3), valImag(3);
1060 fullMatrix<double> vecLeft(3, 3), vecRight(3, 3);
1061 metric.eig(valReal, valImag, vecLeft, vecRight, true);
1062
1063 return std::sqrt(valReal(0) / valReal(2));
1064 }
1065
1066 default:
1067 Msg::Error("wrong dimension for getEigenvaluesMetric function");
1068 return -1;
1069 }
1070 }
1071
pnt(double u,double v,double w,SPoint3 & p) const1072 void MElement::pnt(double u, double v, double w, SPoint3 &p) const
1073 {
1074 double x = 0., y = 0., z = 0.;
1075 double sf[1256];
1076 getShapeFunctions(u, v, w, sf);
1077 for(std::size_t j = 0; j < getNumShapeFunctions(); j++) {
1078 const MVertex *v = getShapeFunctionNode(j);
1079 x += sf[j] * v->x();
1080 y += sf[j] * v->y();
1081 z += sf[j] * v->z();
1082 }
1083 p = SPoint3(x, y, z);
1084 }
1085
pnt(double u,double v,double w,double * p) const1086 void MElement::pnt(double u, double v, double w, double *p) const
1087 {
1088 double x = 0., y = 0., z = 0.;
1089 double sf[1256];
1090 getShapeFunctions(u, v, w, sf);
1091 for(std::size_t j = 0; j < getNumShapeFunctions(); j++) {
1092 const MVertex *v = getShapeFunctionNode(j);
1093 x += sf[j] * v->x();
1094 y += sf[j] * v->y();
1095 z += sf[j] * v->z();
1096 }
1097 p[0] = x;
1098 p[1] = y;
1099 p[2] = z;
1100 }
1101
pnt(const std::vector<double> & sf,SPoint3 & p) const1102 void MElement::pnt(const std::vector<double> &sf, SPoint3 &p) const
1103 {
1104 double x = 0., y = 0., z = 0.;
1105 for(std::size_t j = 0; j < getNumShapeFunctions(); j++) {
1106 const MVertex *v = getShapeFunctionNode(j);
1107 x += sf[j] * v->x();
1108 y += sf[j] * v->y();
1109 z += sf[j] * v->z();
1110 }
1111 p = SPoint3(x, y, z);
1112 }
1113
primaryPnt(double u,double v,double w,SPoint3 & p)1114 void MElement::primaryPnt(double u, double v, double w, SPoint3 &p)
1115 {
1116 double x = 0., y = 0., z = 0.;
1117 double sf[1256];
1118 getShapeFunctions(u, v, w, sf, 1);
1119 for(std::size_t j = 0; j < getNumPrimaryShapeFunctions(); j++) {
1120 const MVertex *v = getShapeFunctionNode(j);
1121 x += sf[j] * v->x();
1122 y += sf[j] * v->y();
1123 z += sf[j] * v->z();
1124 }
1125 p = SPoint3(x, y, z);
1126 }
1127
xyz2uvw(double xyz[3],double uvw[3]) const1128 void MElement::xyz2uvw(double xyz[3], double uvw[3]) const
1129 {
1130 // general Newton routine for the nonlinear case (more efficient
1131 // routines are implemented for simplices, where the basis functions
1132 // are linear)
1133 uvw[0] = uvw[1] = uvw[2] = 0.;
1134
1135 // For high order elements, start from the nearer point
1136 if(getPolynomialOrder() > 2) {
1137 int numNearer = 0;
1138 const MVertex *v = getShapeFunctionNode(0);
1139 double distNearer = (v->x() - xyz[0]) * (v->x() - xyz[0]) +
1140 (v->y() - xyz[1]) * (v->y() - xyz[1]) +
1141 (v->z() - xyz[2]) * (v->z() - xyz[2]);
1142 for(std::size_t i = 1; i < getNumShapeFunctions(); i++) {
1143 const MVertex *v = getShapeFunctionNode(i);
1144 double dist = (v->x() - xyz[0]) * (v->x() - xyz[0]) +
1145 (v->y() - xyz[1]) * (v->y() - xyz[1]) +
1146 (v->z() - xyz[2]) * (v->z() - xyz[2]);
1147 if(dist < distNearer) {
1148 numNearer = i;
1149 distNearer = dist;
1150 }
1151 }
1152 const nodalBasis *nb = getFunctionSpace();
1153 fullMatrix<double> refpnts = nb->getReferenceNodes();
1154 for(int i = 0; i < getDim(); i++) { uvw[i] = refpnts(numNearer, i); }
1155 }
1156
1157 int iter = 1, maxiter = 20;
1158 double error = 1., tol = 1.e-6;
1159
1160 while(error > tol && iter < maxiter) {
1161 double jac[3][3];
1162 if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
1163 double xn = 0., yn = 0., zn = 0.;
1164 double sf[1256];
1165 getShapeFunctions(uvw[0], uvw[1], uvw[2], sf);
1166 for(std::size_t i = 0; i < getNumShapeFunctions(); i++) {
1167 const MVertex *v = getShapeFunctionNode(i);
1168 xn += v->x() * sf[i];
1169 yn += v->y() * sf[i];
1170 zn += v->z() * sf[i];
1171 }
1172 double inv[3][3];
1173 inv3x3(jac, inv);
1174 double un = uvw[0] + inv[0][0] * (xyz[0] - xn) + inv[1][0] * (xyz[1] - yn) +
1175 inv[2][0] * (xyz[2] - zn);
1176 double vn = uvw[1] + inv[0][1] * (xyz[0] - xn) + inv[1][1] * (xyz[1] - yn) +
1177 inv[2][1] * (xyz[2] - zn);
1178 double wn = uvw[2] + inv[0][2] * (xyz[0] - xn) + inv[1][2] * (xyz[1] - yn) +
1179 inv[2][2] * (xyz[2] - zn);
1180 error = sqrt(SQU(un - uvw[0]) + SQU(vn - uvw[1]) + SQU(wn - uvw[2]));
1181 uvw[0] = un;
1182 uvw[1] = vn;
1183 uvw[2] = wn;
1184 iter++;
1185 }
1186 }
1187
movePointFromParentSpaceToElementSpace(double & u,double & v,double & w) const1188 void MElement::movePointFromParentSpaceToElementSpace(double &u, double &v,
1189 double &w) const
1190 {
1191 if(!getParent()) return;
1192 SPoint3 p;
1193 getParent()->pnt(u, v, w, p);
1194 double xyz[3] = {p.x(), p.y(), p.z()};
1195 double uvwE[3];
1196 xyz2uvw(xyz, uvwE);
1197 u = uvwE[0];
1198 v = uvwE[1];
1199 w = uvwE[2];
1200 }
1201
movePointFromElementSpaceToParentSpace(double & u,double & v,double & w) const1202 void MElement::movePointFromElementSpaceToParentSpace(double &u, double &v,
1203 double &w) const
1204 {
1205 if(!getParent()) return;
1206 SPoint3 p;
1207 pnt(u, v, w, p);
1208 double xyz[3] = {p.x(), p.y(), p.z()};
1209 double uvwP[3];
1210 getParent()->xyz2uvw(xyz, uvwP);
1211 u = uvwP[0];
1212 v = uvwP[1];
1213 w = uvwP[2];
1214 }
1215
interpolate(double val[],double u,double v,double w,int stride,int order)1216 double MElement::interpolate(double val[], double u, double v, double w,
1217 int stride, int order)
1218 {
1219 double sum = 0;
1220 int j = 0;
1221 double sf[1256];
1222 getShapeFunctions(u, v, w, sf, order);
1223 for(std::size_t i = 0; i < getNumShapeFunctions(); i++) {
1224 sum += val[j] * sf[i];
1225 j += stride;
1226 }
1227 return sum;
1228 }
1229
interpolateGrad(double val[],double u,double v,double w,double f[],int stride,double invjac[3][3],int order)1230 void MElement::interpolateGrad(double val[], double u, double v, double w,
1231 double f[], int stride, double invjac[3][3],
1232 int order)
1233 {
1234 double dfdu[3] = {0., 0., 0.};
1235 int j = 0;
1236 double gsf[1256][3];
1237 getGradShapeFunctions(u, v, w, gsf, order);
1238 for(std::size_t i = 0; i < getNumShapeFunctions(); i++) {
1239 dfdu[0] += val[j] * gsf[i][0];
1240 dfdu[1] += val[j] * gsf[i][1];
1241 dfdu[2] += val[j] * gsf[i][2];
1242 j += stride;
1243 }
1244 if(invjac) { matvec(invjac, dfdu, f); }
1245 else {
1246 double jac[3][3], inv[3][3];
1247 getJacobian(u, v, w, jac);
1248 inv3x3(jac, inv);
1249 matvec(inv, dfdu, f);
1250 }
1251 }
1252
interpolateCurl(double val[],double u,double v,double w,double f[],int stride,int order)1253 void MElement::interpolateCurl(double val[], double u, double v, double w,
1254 double f[], int stride, int order)
1255 {
1256 double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
1257 getJacobian(u, v, w, jac);
1258 inv3x3(jac, inv);
1259 interpolateGrad(&val[0], u, v, w, fx, stride, inv, order);
1260 interpolateGrad(&val[1], u, v, w, fy, stride, inv, order);
1261 interpolateGrad(&val[2], u, v, w, fz, stride, inv, order);
1262 f[0] = fz[1] - fy[2];
1263 f[1] = -(fz[0] - fx[2]);
1264 f[2] = fy[0] - fx[1];
1265 }
1266
interpolateDiv(double val[],double u,double v,double w,int stride,int order)1267 double MElement::interpolateDiv(double val[], double u, double v, double w,
1268 int stride, int order)
1269 {
1270 double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
1271 getJacobian(u, v, w, jac);
1272 inv3x3(jac, inv);
1273 interpolateGrad(&val[0], u, v, w, fx, stride, inv, order);
1274 interpolateGrad(&val[1], u, v, w, fy, stride, inv, order);
1275 interpolateGrad(&val[2], u, v, w, fz, stride, inv, order);
1276 return fx[0] + fy[1] + fz[2];
1277 }
1278
integrate(double val[],int pOrder,int stride,int order)1279 double MElement::integrate(double val[], int pOrder, int stride, int order)
1280 {
1281 int npts;
1282 IntPt *gp;
1283 getIntegrationPoints(pOrder, &npts, &gp);
1284 double sum = 0;
1285 for(int i = 0; i < npts; i++) {
1286 double u = gp[i].pt[0];
1287 double v = gp[i].pt[1];
1288 double w = gp[i].pt[2];
1289 double weight = gp[i].weight;
1290 double detuvw = getJacobianDeterminant(u, v, w);
1291 sum += interpolate(val, u, v, w, stride, order) * weight * detuvw;
1292 }
1293 return sum;
1294 }
1295
integrateCirc(double val[],int edge,int pOrder,int order)1296 double MElement::integrateCirc(double val[], int edge, int pOrder, int order)
1297 {
1298 if(edge > getNumEdges() - 1) {
1299 Msg::Error("No edge %d for this element", edge);
1300 return 0;
1301 }
1302
1303 std::vector<MVertex *> v;
1304 getEdgeVertices(edge, v);
1305 MElementFactory f;
1306 int type = ElementType::getType(TYPE_LIN, getPolynomialOrder());
1307 MElement *ee = f.create(type, v);
1308
1309 double intv[3];
1310 for(int i = 0; i < 3; i++) {
1311 intv[i] = ee->integrate(&val[i], pOrder, 3, order);
1312 }
1313 delete ee;
1314
1315 double t[3] = {v[1]->x() - v[0]->x(), v[1]->y() - v[0]->y(),
1316 v[1]->z() - v[0]->z()};
1317 norme(t);
1318
1319 return prosca(t, intv);
1320 }
1321
integrateFlux(double val[],int face,int pOrder,int order)1322 double MElement::integrateFlux(double val[], int face, int pOrder, int order)
1323 {
1324 if(face > getNumFaces() - 1) {
1325 Msg::Error("No face %d for this element", face);
1326 return 0;
1327 }
1328 std::vector<MVertex *> v;
1329 getFaceVertices(face, v);
1330 MElementFactory f;
1331 int type = 0;
1332 switch(getType()) {
1333 case TYPE_TRI:
1334 case TYPE_TET:
1335 case TYPE_QUA:
1336 case TYPE_HEX:
1337 type = ElementType::getType(getType(), getPolynomialOrder());
1338 break;
1339 case TYPE_PYR:
1340 if(face < 4)
1341 type = ElementType::getType(TYPE_TRI, getPolynomialOrder());
1342 else
1343 type = ElementType::getType(TYPE_QUA, getPolynomialOrder());
1344 break;
1345 case TYPE_PRI:
1346 if(face < 2)
1347 type = ElementType::getType(TYPE_TRI, getPolynomialOrder());
1348 else
1349 type = ElementType::getType(TYPE_QUA, getPolynomialOrder());
1350 break;
1351 default: type = 0; break;
1352 }
1353 MElement *fe = f.create(type, v);
1354
1355 double intv[3];
1356 for(int i = 0; i < 3; i++) {
1357 intv[i] = fe->integrate(&val[i], pOrder, 3, order);
1358 }
1359 delete fe;
1360
1361 double n[3];
1362 normal3points(v[0]->x(), v[0]->y(), v[0]->z(), v[1]->x(), v[1]->y(),
1363 v[1]->z(), v[2]->x(), v[2]->y(), v[2]->z(), n);
1364
1365 return prosca(n, intv);
1366 }
1367
writeMSH2(FILE * fp,double version,bool binary,int num,int elementary,int physical,int parentNum,int dom1Num,int dom2Num,std::vector<short> * ghosts)1368 void MElement::writeMSH2(FILE *fp, double version, bool binary, int num,
1369 int elementary, int physical, int parentNum,
1370 int dom1Num, int dom2Num, std::vector<short> *ghosts)
1371 {
1372 int type = getTypeForMSH();
1373
1374 if(!type) return;
1375
1376 int n = getNumVerticesForMSH();
1377 int par = (parentNum) ? 1 : 0;
1378 int dom = (dom1Num) ? 2 : 0;
1379 bool poly = (type == MSH_POLYG_ || type == MSH_POLYH_ || type == MSH_POLYG_B);
1380
1381 // if polygon loop over children (triangles and tets)
1382 if(CTX::instance()->mesh.saveTri) {
1383 if(poly) {
1384 for(int i = 0; i < getNumChildren(); i++) {
1385 MElement *t = getChild(i);
1386 t->writeMSH2(fp, version, binary, num++, elementary, physical, 0, 0, 0,
1387 ghosts);
1388 }
1389 return;
1390 }
1391 if(type == MSH_TRI_B) {
1392 MTriangle *t = new MTriangle(getVertex(0), getVertex(1), getVertex(2));
1393 t->writeMSH2(fp, version, binary, num++, elementary, physical, 0, 0, 0,
1394 ghosts);
1395 delete t;
1396 return;
1397 }
1398 if(type == MSH_LIN_B || type == MSH_LIN_C) {
1399 MLine *l = new MLine(getVertex(0), getVertex(1));
1400 l->writeMSH2(fp, version, binary, num++, elementary, physical, 0, 0, 0,
1401 ghosts);
1402 delete l;
1403 return;
1404 }
1405 }
1406
1407 if(CTX::instance()->mesh.preserveNumberingMsh2) num = (int)_num;
1408
1409 if(!binary) {
1410 fprintf(fp, "%d %d", num ? num : (int)_num, type);
1411 if(version < 2.0)
1412 fprintf(fp, " %d %d %d", abs(physical), elementary, n);
1413 else if(version < 2.2)
1414 fprintf(fp, " %d %d %d", abs(physical), elementary, _partition);
1415 else if(!_partition && !par && !dom)
1416 fprintf(fp, " %d %d %d", 2 + par + dom, abs(physical), elementary);
1417 else if(!ghosts)
1418 fprintf(fp, " %d %d %d 1 %d", 4 + par + dom, abs(physical), elementary,
1419 _partition);
1420 else {
1421 int numGhosts = ghosts->size();
1422 fprintf(fp, " %d %d %d %d %d", 4 + numGhosts + par + dom, abs(physical),
1423 elementary, 1 + numGhosts, _partition);
1424 for(std::size_t i = 0; i < ghosts->size(); i++)
1425 fprintf(fp, " %d", -(*ghosts)[i]);
1426 }
1427 if(version >= 2.0 && par) fprintf(fp, " %d", parentNum);
1428 if(version >= 2.0 && dom) fprintf(fp, " %d %d", dom1Num, dom2Num);
1429 if(version >= 2.0 && poly) fprintf(fp, " %d", n);
1430 }
1431 else {
1432 int numTags, numGhosts = 0;
1433 if(!_partition)
1434 numTags = 2;
1435 else if(!ghosts)
1436 numTags = 4;
1437 else {
1438 numGhosts = ghosts->size();
1439 numTags = 4 + numGhosts;
1440 }
1441 numTags += par;
1442 // we write elements in blobs of single elements; this will lead
1443 // to suboptimal reads, but it's much simpler when the number of
1444 // tags change from element to element (third-party codes can
1445 // still write MSH file optimized for reading speed, by grouping
1446 // elements with the same number of tags in blobs)
1447 int blob[60] = {
1448 type, 1, numTags, num ? num : (int)_num,
1449 abs(physical), elementary, 1 + numGhosts, _partition};
1450 if(ghosts)
1451 for(int i = 0; i < numGhosts; i++) blob[8 + i] = -(*ghosts)[i];
1452 if(par) blob[8 + numGhosts] = parentNum;
1453 if(poly) Msg::Error("Unable to write polygons/polyhedra in binary files.");
1454 fwrite(blob, sizeof(int), 4 + numTags, fp);
1455 }
1456
1457 if(physical < 0) reverse();
1458
1459 std::vector<int> verts;
1460 getVerticesIdForMSH(verts);
1461
1462 if(!binary) {
1463 for(int i = 0; i < n; i++) fprintf(fp, " %d", verts[i]);
1464 fprintf(fp, "\n");
1465 }
1466 else {
1467 fwrite(&verts[0], sizeof(int), n, fp);
1468 }
1469
1470 if(physical < 0) reverse();
1471 }
1472
writeMSH3(FILE * fp,bool binary,int entity,std::vector<short> * ghosts)1473 void MElement::writeMSH3(FILE *fp, bool binary, int entity,
1474 std::vector<short> *ghosts)
1475 {
1476 int num = (int)getNum();
1477 int type = getTypeForMSH();
1478 if(!type) return;
1479
1480 std::vector<int> verts;
1481 getVerticesIdForMSH(verts);
1482
1483 // FIXME: once we create elements using their own interpretion of data, we
1484 // should move this also into each element base class
1485 std::vector<int> data;
1486 data.insert(data.end(), verts.begin(), verts.end());
1487 if(getParent()) data.push_back(getParent()->getNum());
1488 if(getPartition()) {
1489 if(ghosts) {
1490 data.push_back(1 + ghosts->size());
1491 data.push_back(getPartition());
1492 data.insert(data.end(), ghosts->begin(), ghosts->end());
1493 }
1494 else {
1495 data.push_back(1);
1496 data.push_back(getPartition());
1497 }
1498 }
1499 int numData = data.size();
1500
1501 if(!binary) {
1502 fprintf(fp, "%d %d %d %d", num, type, entity, numData);
1503 for(int i = 0; i < numData; i++) fprintf(fp, " %d", data[i]);
1504 fprintf(fp, "\n");
1505 }
1506 else {
1507 fwrite(&num, sizeof(int), 1, fp);
1508 fwrite(&type, sizeof(int), 1, fp);
1509 fwrite(&entity, sizeof(int), 1, fp);
1510 fwrite(&numData, sizeof(int), 1, fp);
1511 fwrite(&data[0], sizeof(int), numData, fp);
1512 }
1513 }
1514
writePOS(FILE * fp,bool printElementary,bool printElementNumber,bool printSICN,bool printSIGE,bool printGamma,bool printDisto,double scalingFactor,int elementary)1515 void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
1516 bool printSICN, bool printSIGE, bool printGamma,
1517 bool printDisto, double scalingFactor, int elementary)
1518 {
1519 const char *str = getStringForPOS();
1520 if(!str) return;
1521
1522 int n = getNumVertices();
1523 fprintf(fp, "%s(", str);
1524 for(int i = 0; i < n; i++) {
1525 if(i) fprintf(fp, ",");
1526 fprintf(fp, "%g,%g,%g", getVertex(i)->x() * scalingFactor,
1527 getVertex(i)->y() * scalingFactor,
1528 getVertex(i)->z() * scalingFactor);
1529 }
1530 fprintf(fp, "){");
1531 bool first = true;
1532 if(printElementary) {
1533 for(int i = 0; i < n; i++) {
1534 if(first)
1535 first = false;
1536 else
1537 fprintf(fp, ",");
1538 fprintf(fp, "%d", elementary);
1539 }
1540 }
1541 if(printElementNumber) {
1542 for(int i = 0; i < n; i++) {
1543 if(first)
1544 first = false;
1545 else
1546 fprintf(fp, ",");
1547 fprintf(fp, "%lu", getNum());
1548 }
1549 }
1550 if(printSICN) {
1551 double sICNMin = minSICNShapeMeasure();
1552 for(int i = 0; i < n; i++) {
1553 if(first)
1554 first = false;
1555 else
1556 fprintf(fp, ",");
1557 fprintf(fp, "%g", sICNMin);
1558 }
1559 }
1560 if(printSIGE) {
1561 double sIGEMin = minSIGEShapeMeasure();
1562 for(int i = 0; i < n; i++) {
1563 if(first)
1564 first = false;
1565 else
1566 fprintf(fp, ",");
1567 fprintf(fp, "%g", sIGEMin);
1568 }
1569 }
1570 if(printGamma) {
1571 double gamma = gammaShapeMeasure();
1572 for(int i = 0; i < n; i++) {
1573 if(first)
1574 first = false;
1575 else
1576 fprintf(fp, ",");
1577 fprintf(fp, "%g", gamma);
1578 }
1579 }
1580 if(printDisto) {
1581 double disto = distoShapeMeasure();
1582 for(int i = 0; i < n; i++) {
1583 if(first)
1584 first = false;
1585 else
1586 fprintf(fp, ",");
1587 fprintf(fp, "%g", disto);
1588 }
1589 }
1590 fprintf(fp, "};\n");
1591 }
1592
writeSTL(FILE * fp,bool binary,double scalingFactor)1593 void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
1594 {
1595 if(getType() != TYPE_TRI && getType() != TYPE_QUA) return;
1596 int qid[3] = {0, 2, 3};
1597 SVector3 n = getFace(0).normal();
1598 if(!binary) {
1599 fprintf(fp, "facet normal %g %g %g\n", n[0], n[1], n[2]);
1600 fprintf(fp, " outer loop\n");
1601 for(int j = 0; j < 3; j++)
1602 fprintf(
1603 fp, " vertex %.16g %.16g %.16g\n", getVertex(j)->x() * scalingFactor,
1604 getVertex(j)->y() * scalingFactor, getVertex(j)->z() * scalingFactor);
1605 fprintf(fp, " endloop\n");
1606 fprintf(fp, "endfacet\n");
1607 if(getNumVertices() == 4) {
1608 fprintf(fp, "facet normal %g %g %g\n", n[0], n[1], n[2]);
1609 fprintf(fp, " outer loop\n");
1610 for(int j = 0; j < 3; j++)
1611 fprintf(fp, " vertex %.16g %.16g %.16g\n",
1612 getVertex(qid[j])->x() * scalingFactor,
1613 getVertex(qid[j])->y() * scalingFactor,
1614 getVertex(qid[j])->z() * scalingFactor);
1615 fprintf(fp, " endloop\n");
1616 fprintf(fp, "endfacet\n");
1617 }
1618 }
1619 else {
1620 char data[50];
1621 float *coords = (float *)data;
1622 coords[0] = (float)n[0];
1623 coords[1] = (float)n[1];
1624 coords[2] = (float)n[2];
1625 for(int j = 0; j < 3; j++) {
1626 coords[3 + 3 * j] = (float)(getVertex(j)->x() * scalingFactor);
1627 coords[3 + 3 * j + 1] = (float)(getVertex(j)->y() * scalingFactor);
1628 coords[3 + 3 * j + 2] = (float)(getVertex(j)->z() * scalingFactor);
1629 }
1630 data[48] = data[49] = 0;
1631 fwrite(data, sizeof(char), 50, fp);
1632 if(getNumVertices() == 4) {
1633 for(int j = 0; j < 3; j++) {
1634 coords[3 + 3 * j] = (float)(getVertex(qid[j])->x() * scalingFactor);
1635 coords[3 + 3 * j + 1] = (float)(getVertex(qid[j])->y() * scalingFactor);
1636 coords[3 + 3 * j + 2] = (float)(getVertex(qid[j])->z() * scalingFactor);
1637 }
1638 fwrite(data, sizeof(char), 50, fp);
1639 }
1640 }
1641 }
1642
writeX3D(FILE * fp,double scalingFactor)1643 void MElement::writeX3D(FILE *fp, double scalingFactor)
1644 {
1645 if(getType() != TYPE_TRI && getType() != TYPE_QUA) return;
1646 int qid[3] = {0, 2, 3};
1647
1648 for(int j = 0; j < 3; j++)
1649 fprintf(fp, "%g %g %g\n", getVertex(j)->x() * scalingFactor,
1650 getVertex(j)->y() * scalingFactor,
1651 getVertex(j)->z() * scalingFactor);
1652 if(getNumVertices() == 4) {
1653 for(int j = 0; j < 3; j++) {
1654 fprintf(fp, "%g %g %g\n", getVertex(qid[j])->x() * scalingFactor,
1655 getVertex(qid[j])->y() * scalingFactor,
1656 getVertex(qid[j])->z() * scalingFactor);
1657 }
1658 }
1659 }
1660
writePLY2(FILE * fp)1661 void MElement::writePLY2(FILE *fp)
1662 {
1663 fprintf(fp, "3 ");
1664 for(std::size_t i = 0; i < getNumVertices(); i++)
1665 fprintf(fp, " %ld", getVertex(i)->getIndex() - 1);
1666 fprintf(fp, "\n");
1667 }
1668
writeVRML(FILE * fp)1669 void MElement::writeVRML(FILE *fp)
1670 {
1671 for(std::size_t i = 0; i < getNumVertices(); i++)
1672 fprintf(fp, "%ld,", getVertex(i)->getIndex() - 1);
1673 fprintf(fp, "-1,\n");
1674 }
1675
writeTOCHNOG(FILE * fp,int num)1676 void MElement::writeTOCHNOG(FILE *fp, int num)
1677 {
1678 const char *str = getStringForTOCHNOG();
1679 if(!str) return;
1680
1681 int n = getNumVertices();
1682 fprintf(fp, "element %d %s ", num, str);
1683 for(int i = 0; i < n; i++) {
1684 fprintf(fp, " %ld", getVertexTOCHNOG(i)->getIndex());
1685 }
1686 fprintf(fp, "\n");
1687 }
1688
writeVTK(FILE * fp,bool binary,bool bigEndian)1689 void MElement::writeVTK(FILE *fp, bool binary, bool bigEndian)
1690 {
1691 if(!getTypeForVTK()) return;
1692
1693 int n = getNumVertices();
1694 if(binary) {
1695 int verts[60];
1696 verts[0] = n;
1697 for(int i = 0; i < n; i++)
1698 verts[i + 1] = (int)getVertexVTK(i)->getIndex() - 1;
1699 // VTK always expects big endian binary data
1700 if(!bigEndian) SwapBytes((char *)verts, sizeof(int), n + 1);
1701 fwrite(verts, sizeof(int), n + 1, fp);
1702 }
1703 else {
1704 fprintf(fp, "%d", n);
1705 for(int i = 0; i < n; i++)
1706 fprintf(fp, " %ld", getVertexVTK(i)->getIndex() - 1);
1707 fprintf(fp, "\n");
1708 }
1709 }
1710
writeMATLAB(FILE * fp,int filetype,int elementary,int physical,bool binary)1711 void MElement::writeMATLAB(FILE *fp, int filetype, int elementary, int physical,
1712 bool binary)
1713 {
1714 // Matlab use the same names as MSH
1715 if(!getTypeForMSH()) return;
1716 if(binary) {
1717 Msg::Warning(
1718 "Binary format not available for Matlab, saving into ASCII format");
1719 binary = false;
1720 }
1721
1722 // Simple version
1723 if(filetype == 0) {
1724 int n = getNumVertices();
1725 for(int i = 0; i < n; i++)
1726 fprintf(fp, " %ld", getVertexMATLAB(i)->getIndex());
1727 fprintf(fp, ";\n");
1728 }
1729 // same as load_gmsh2.m
1730 if(filetype == 1) {
1731 if(physical < 0) reverse();
1732
1733 for(std::size_t i = 0; i < getNumVertices(); i++)
1734 fprintf(fp, " %ld", getVertex(i)->getIndex());
1735 fprintf(fp, " %d\n", physical ? abs(physical) : elementary);
1736
1737 if(physical < 0) reverse();
1738 }
1739 }
1740
writeUNV(FILE * fp,int num,int elementary,int physical)1741 void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
1742 {
1743 int type = getTypeForUNV();
1744 if(!type) return;
1745
1746 int n = getNumVertices();
1747 int physical_property = elementary;
1748 int material_property = abs(physical);
1749 int color = 7;
1750 fprintf(fp, "%10d%10d%10d%10d%10d%10d\n", num ? num : (int)_num, type,
1751 physical_property, material_property, color, n);
1752 if(type == 21 || type == 24) // linear beam or parabolic beam
1753 fprintf(fp, "%10d%10d%10d\n", 0, 0, 0);
1754
1755 if(physical < 0) reverse();
1756
1757 for(int k = 0; k < n; k++) {
1758 fprintf(fp, "%10ld", getVertexUNV(k)->getIndex());
1759 if(k % 8 == 7) fprintf(fp, "\n");
1760 }
1761 if(n - 1 % 8 != 7) fprintf(fp, "\n");
1762
1763 if(physical < 0) reverse();
1764 }
1765
writeMESH(FILE * fp,int elementTagType,int elementary,int physical)1766 void MElement::writeMESH(FILE *fp, int elementTagType, int elementary,
1767 int physical)
1768 {
1769 if(physical < 0) reverse();
1770
1771 for(std::size_t i = 0; i < getNumVertices(); i++)
1772 if(getTypeForMSH() == MSH_TET_10 && i == 8)
1773 fprintf(fp, " %ld", getVertex(9)->getIndex());
1774 else if(getTypeForMSH() == MSH_TET_10 && i == 9)
1775 fprintf(fp, " %ld", getVertex(8)->getIndex());
1776 else
1777 fprintf(fp, " %ld", getVertex(i)->getIndex());
1778 fprintf(fp, " %d\n",
1779 (elementTagType == 3) ? _partition :
1780 (elementTagType == 2) ? abs(physical) :
1781 elementary);
1782
1783 if(physical < 0) reverse();
1784 }
1785
writeNEU(FILE * fp,unsigned gambitType,int idAdjust,int phys)1786 void MElement::writeNEU(FILE *fp, unsigned gambitType, int idAdjust, int phys)
1787 {
1788 if(phys < 0) reverse();
1789
1790 fprintf(fp, "%8lu %2d %2lu ", _num - idAdjust, gambitType, getNumVertices());
1791 for(std::size_t i = 0; i < getNumVertices(); ++i) {
1792 if(i == 7) fprintf(fp, "\n ");
1793 fprintf(fp, "%8ld", getVertexNEU(i)->getIndex());
1794 }
1795 fprintf(fp, "\n");
1796
1797 if(phys < 0) reverse();
1798 }
1799
writeIR3(FILE * fp,int elementTagType,int num,int elementary,int physical)1800 void MElement::writeIR3(FILE *fp, int elementTagType, int num, int elementary,
1801 int physical)
1802 {
1803 if(physical < 0) reverse();
1804
1805 int numVert = getNumVertices();
1806 fprintf(fp, "%d %d %d", num,
1807 (elementTagType == 3) ? _partition :
1808 (elementTagType == 2) ? abs(physical) :
1809 elementary,
1810 numVert);
1811 for(int i = 0; i < numVert; i++)
1812 fprintf(fp, " %ld", getVertex(i)->getIndex());
1813 fprintf(fp, "\n");
1814
1815 if(physical < 0) reverse();
1816 }
1817
writeBDF(FILE * fp,int format,int elementTagType,int elementary,int physical)1818 void MElement::writeBDF(FILE *fp, int format, int elementTagType,
1819 int elementary, int physical)
1820 {
1821 const char *str = getStringForBDF();
1822 if(!str) return;
1823
1824 int n = getNumVertices();
1825 const char *cont[4] = {"E", "F", "G", "H"};
1826 int ncont = 0;
1827
1828 if(physical < 0) reverse();
1829
1830 int tag = (elementTagType == 3) ? _partition :
1831 (elementTagType == 2) ? abs(physical) :
1832 elementary;
1833
1834 if(format == 0) { // free field format
1835 fprintf(fp, "%s,%lu,%d", str, _num, tag);
1836 for(int i = 0; i < n; i++) {
1837 fprintf(fp, ",%ld", getVertexBDF(i)->getIndex());
1838 if(i != n - 1 && !((i + 3) % 8)) {
1839 fprintf(fp, ",+%s%lu\n+%s%lu", cont[ncont], _num, cont[ncont], _num);
1840 ncont++;
1841 }
1842 }
1843 if(n == 2) // CBAR
1844 fprintf(fp, ",0.,0.,0.");
1845 fprintf(fp, "\n");
1846 }
1847 else if(format == 1) { // small field format
1848 fprintf(fp, "%-8s%-8lu%-8d", str, _num, tag);
1849 for(int i = 0; i < n; i++) {
1850 fprintf(fp, "%-8ld", getVertexBDF(i)->getIndex());
1851 if(i != n - 1 && !((i + 3) % 8)) {
1852 fprintf(fp, "+%s%-6lu\n+%s%-6lu", cont[ncont], _num, cont[ncont], _num);
1853 ncont++;
1854 }
1855 }
1856 if(n == 2) // CBAR
1857 fprintf(fp, "%-8s%-8s%-8s", "0.", "0.", "0.");
1858 fprintf(fp, "\n");
1859 }
1860 else { // large field format
1861 fprintf(fp, "%-8s%-8lu%-8d", str, _num, tag);
1862 for(int i = 0; i < n; i++) {
1863 fprintf(fp, "%-8ld", getVertexBDF(i)->getIndex());
1864 if(i != n - 1 && !((i + 3) % 8)) {
1865 fprintf(fp, "\n ");
1866 ncont++;
1867 }
1868 }
1869 if(n == 2) // CBAR
1870 fprintf(fp, "%-8s%-8s%-8s", "0.", "0.", "0.");
1871 fprintf(fp, "\n");
1872 }
1873
1874 if(physical < 0) reverse();
1875 }
1876
writeDIFF(FILE * fp,int num,bool binary,int physical)1877 void MElement::writeDIFF(FILE *fp, int num, bool binary, int physical)
1878 {
1879 const char *str = getStringForDIFF();
1880 if(!str) return;
1881
1882 if(physical < 0) reverse();
1883
1884 int n = getNumVertices();
1885 if(binary) {
1886 // TODO
1887 }
1888 else {
1889 fprintf(fp, "%d %s %d ", num, str, abs(physical));
1890 for(int i = 0; i < n; i++)
1891 fprintf(fp, " %ld", getVertexDIFF(i)->getIndex());
1892 fprintf(fp, "\n");
1893 }
1894
1895 if(physical < 0) reverse();
1896 }
1897
writeINP(FILE * fp,int num)1898 void MElement::writeINP(FILE *fp, int num)
1899 {
1900 fprintf(fp, "%d, ", num);
1901 int n = getNumVertices();
1902 for(int i = 0; i < n; i++) {
1903 fprintf(fp, "%ld", getVertexINP(i)->getIndex());
1904 if(i != n - 1) {
1905 fprintf(fp, ", ");
1906 if(i && !((i + 2) % 16)) fprintf(fp, "\n");
1907 }
1908 }
1909 fprintf(fp, "\n");
1910 }
1911
writeKEY(FILE * fp,int pid,int num)1912 void MElement::writeKEY(FILE *fp, int pid, int num)
1913 {
1914 fprintf(fp, "%d, %d, ", num, pid);
1915 int n = getNumVertices();
1916 int nid[64];
1917 int i;
1918 for(i = 0; i < n; i++) nid[i] = (int)getVertexKEY(i)->getIndex();
1919 if(getDim() == 3) {
1920 if(n == 4) { /* tet4, repeating n4 */
1921 nid[7] = nid[6] = nid[5] = nid[4] = nid[3];
1922 }
1923 else if(n == 6) { /* penta6, n8=n7 & n4=n3 */
1924 nid[7] = nid[6] = nid[5];
1925 nid[5] = nid[4];
1926 }
1927 if(n < 8) n = 8;
1928 }
1929 else if(getDim() == 2) {
1930 if(n == 3) { /* 3-node shell */
1931 nid[3] = nid[2];
1932 n++;
1933 }
1934 else if(n == 6) { /* 6-node shell */
1935 nid[7] = nid[6] = nid[5];
1936 nid[5] = nid[4];
1937 nid[4] = nid[3];
1938 nid[3] = nid[2];
1939 n = 8;
1940 }
1941 }
1942 else if(getDim() == 1) {
1943 if(n == 3) { /* elbow, write the third node on the next line */
1944 nid[8] = nid[2];
1945 for(i = 2; i < 8; i++) nid[i] = 0;
1946 n = 9;
1947 }
1948 }
1949 for(i = 0; i < n; i++) {
1950 fprintf(fp, "%d", nid[i]);
1951 if(i != n - 1) {
1952 fprintf(fp, ", ");
1953 if(!((i + 2) % 10)) fprintf(fp, "\n");
1954 }
1955 }
1956 fprintf(fp, "\n");
1957 }
1958
writeSU2(FILE * fp,int num)1959 void MElement::writeSU2(FILE *fp, int num)
1960 {
1961 fprintf(fp, "%d ", getTypeForVTK());
1962 for(std::size_t i = 0; i < getNumVertices(); i++)
1963 fprintf(fp, "%ld ", getVertexVTK(i)->getIndex() - 1);
1964 if(num >= 0)
1965 fprintf(fp, "%d\n", num);
1966 else
1967 fprintf(fp, "\n");
1968 }
1969
getInfoMSH(const int typeMSH,const char ** const name)1970 unsigned int MElement::getInfoMSH(const int typeMSH, const char **const name)
1971 {
1972 switch(typeMSH) {
1973 case MSH_PNT:
1974 if(name) *name = "Point";
1975 return 1;
1976 case MSH_LIN_1:
1977 if(name) *name = "Line 1";
1978 return 1;
1979 case MSH_LIN_2:
1980 if(name) *name = "Line 2";
1981 return 2;
1982 case MSH_LIN_3:
1983 if(name) *name = "Line 3";
1984 return 2 + 1;
1985 case MSH_LIN_4:
1986 if(name) *name = "Line 4";
1987 return 2 + 2;
1988 case MSH_LIN_5:
1989 if(name) *name = "Line 5";
1990 return 2 + 3;
1991 case MSH_LIN_6:
1992 if(name) *name = "Line 6";
1993 return 2 + 4;
1994 case MSH_LIN_7:
1995 if(name) *name = "Line 7";
1996 return 2 + 5;
1997 case MSH_LIN_8:
1998 if(name) *name = "Line 8";
1999 return 2 + 6;
2000 case MSH_LIN_9:
2001 if(name) *name = "Line 9";
2002 return 2 + 7;
2003 case MSH_LIN_10:
2004 if(name) *name = "Line 10";
2005 return 2 + 8;
2006 case MSH_LIN_11:
2007 if(name) *name = "Line 11";
2008 return 2 + 9;
2009 case MSH_LIN_B:
2010 if(name) *name = "Line Border";
2011 return 2;
2012 case MSH_LIN_C:
2013 if(name) *name = "Line Child";
2014 return 2;
2015 case MSH_TRI_1:
2016 if(name) *name = "Triangle 1";
2017 return 1;
2018 case MSH_TRI_3:
2019 if(name) *name = "Triangle 3";
2020 return 3;
2021 case MSH_TRI_6:
2022 if(name) *name = "Triangle 6";
2023 return 3 + 3;
2024 case MSH_TRI_9:
2025 if(name) *name = "Triangle 9";
2026 return 3 + 6;
2027 case MSH_TRI_10:
2028 if(name) *name = "Triangle 10";
2029 return 3 + 6 + 1;
2030 case MSH_TRI_12:
2031 if(name) *name = "Triangle 12";
2032 return 3 + 9;
2033 case MSH_TRI_15:
2034 if(name) *name = "Triangle 15";
2035 return 3 + 9 + 3;
2036 case MSH_TRI_15I:
2037 if(name) *name = "Triangle 15I";
2038 return 3 + 12;
2039 case MSH_TRI_21:
2040 if(name) *name = "Triangle 21";
2041 return 3 + 12 + 6;
2042 case MSH_TRI_28:
2043 if(name) *name = "Triangle 28";
2044 return 3 + 15 + 10;
2045 case MSH_TRI_36:
2046 if(name) *name = "Triangle 36";
2047 return 3 + 18 + 15;
2048 case MSH_TRI_45:
2049 if(name) *name = "Triangle 45";
2050 return 3 + 21 + 21;
2051 case MSH_TRI_55:
2052 if(name) *name = "Triangle 55";
2053 return 3 + 24 + 28;
2054 case MSH_TRI_66:
2055 if(name) *name = "Triangle 66";
2056 return 3 + 27 + 36;
2057 case MSH_TRI_18:
2058 if(name) *name = "Triangle 18";
2059 return 3 + 15;
2060 case MSH_TRI_21I:
2061 if(name) *name = "Triangle 21I";
2062 return 3 + 18;
2063 case MSH_TRI_24:
2064 if(name) *name = "Triangle 24";
2065 return 3 + 21;
2066 case MSH_TRI_27:
2067 if(name) *name = "Triangle 27";
2068 return 3 + 24;
2069 case MSH_TRI_30:
2070 if(name) *name = "Triangle 30";
2071 return 3 + 27;
2072 case MSH_TRI_B:
2073 if(name) *name = "Triangle Border";
2074 return 3;
2075 case MSH_QUA_1:
2076 if(name) *name = "Quadrilateral 1";
2077 return 1;
2078 case MSH_QUA_4:
2079 if(name) *name = "Quadrilateral 4";
2080 return 4;
2081 case MSH_QUA_8:
2082 if(name) *name = "Quadrilateral 8";
2083 return 4 + 4;
2084 case MSH_QUA_9:
2085 if(name) *name = "Quadrilateral 9";
2086 return 9;
2087 case MSH_QUA_16:
2088 if(name) *name = "Quadrilateral 16";
2089 return 16;
2090 case MSH_QUA_25:
2091 if(name) *name = "Quadrilateral 25";
2092 return 25;
2093 case MSH_QUA_36:
2094 if(name) *name = "Quadrilateral 36";
2095 return 36;
2096 case MSH_QUA_49:
2097 if(name) *name = "Quadrilateral 49";
2098 return 49;
2099 case MSH_QUA_64:
2100 if(name) *name = "Quadrilateral 64";
2101 return 64;
2102 case MSH_QUA_81:
2103 if(name) *name = "Quadrilateral 81";
2104 return 81;
2105 case MSH_QUA_100:
2106 if(name) *name = "Quadrilateral 100";
2107 return 100;
2108 case MSH_QUA_121:
2109 if(name) *name = "Quadrilateral 121";
2110 return 121;
2111 case MSH_QUA_12:
2112 if(name) *name = "Quadrilateral 12";
2113 return 12;
2114 case MSH_QUA_16I:
2115 if(name) *name = "Quadrilateral 16I";
2116 return 16;
2117 case MSH_QUA_20:
2118 if(name) *name = "Quadrilateral 20";
2119 return 20;
2120 case MSH_QUA_24:
2121 if(name) *name = "Quadrilateral 24";
2122 return 24;
2123 case MSH_QUA_28:
2124 if(name) *name = "Quadrilateral 28";
2125 return 28;
2126 case MSH_QUA_32:
2127 if(name) *name = "Quadrilateral 32";
2128 return 32;
2129 case MSH_QUA_36I:
2130 if(name) *name = "Quadrilateral 36I";
2131 return 36;
2132 case MSH_QUA_40:
2133 if(name) *name = "Quadrilateral 40";
2134 return 40;
2135 case MSH_POLYG_:
2136 if(name) *name = "Polygon";
2137 return 0;
2138 case MSH_POLYG_B:
2139 if(name) *name = "Polygon Border";
2140 return 0;
2141 case MSH_TET_1:
2142 if(name) *name = "Tetrahedron 1";
2143 return 1;
2144 case MSH_TET_4:
2145 if(name) *name = "Tetrahedron 4";
2146 return 4;
2147 case MSH_TET_10:
2148 if(name) *name = "Tetrahedron 10";
2149 return 4 + 6;
2150 case MSH_TET_20:
2151 if(name) *name = "Tetrahedron 20";
2152 return 4 + 12 + 4;
2153 case MSH_TET_35:
2154 if(name) *name = "Tetrahedron 35";
2155 return 4 + 18 + 12 + 1;
2156 case MSH_TET_56:
2157 if(name) *name = "Tetrahedron 56";
2158 return 4 + 24 + 24 + 4;
2159 case MSH_TET_84:
2160 if(name) *name = "Tetrahedron 84";
2161 return (7 * 8 * 9) / 6;
2162 case MSH_TET_120:
2163 if(name) *name = "Tetrahedron 120";
2164 return (8 * 9 * 10) / 6;
2165 case MSH_TET_165:
2166 if(name) *name = "Tetrahedron 165";
2167 return (9 * 10 * 11) / 6;
2168 case MSH_TET_220:
2169 if(name) *name = "Tetrahedron 220";
2170 return (10 * 11 * 12) / 6;
2171 case MSH_TET_286:
2172 if(name) *name = "Tetrahedron 286";
2173 return (11 * 12 * 13) / 6;
2174 case MSH_TET_16:
2175 if(name) *name = "Tetrahedron 16";
2176 return 4 + 6 * 2;
2177 case MSH_TET_22:
2178 if(name) *name = "Tetrahedron 22";
2179 return 4 + 6 * 3;
2180 case MSH_TET_28:
2181 if(name) *name = "Tetrahedron 28";
2182 return 4 + 6 * 4;
2183 case MSH_TET_34:
2184 if(name) *name = "Tetrahedron 34";
2185 return 4 + 6 * 5;
2186 case MSH_TET_40:
2187 if(name) *name = "Tetrahedron 40";
2188 return 4 + 6 * 6;
2189 case MSH_TET_46:
2190 if(name) *name = "Tetrahedron 46";
2191 return 4 + 6 * 7;
2192 case MSH_TET_52:
2193 if(name) *name = "Tetrahedron 52";
2194 return 4 + 6 * 8;
2195 case MSH_TET_58:
2196 if(name) *name = "Tetrahedron 58";
2197 return 4 + 6 * 9;
2198 case MSH_HEX_1:
2199 if(name) *name = "Hexahedron 1";
2200 return 1;
2201 case MSH_HEX_8:
2202 if(name) *name = "Hexahedron 8";
2203 return 8;
2204 case MSH_HEX_20:
2205 if(name) *name = "Hexahedron 20";
2206 return 8 + 12;
2207 case MSH_HEX_27:
2208 if(name) *name = "Hexahedron 27";
2209 return 8 + 12 + 6 + 1;
2210 case MSH_HEX_64:
2211 if(name) *name = "Hexahedron 64";
2212 return 64;
2213 case MSH_HEX_125:
2214 if(name) *name = "Hexahedron 125";
2215 return 125;
2216 case MSH_HEX_216:
2217 if(name) *name = "Hexahedron 216";
2218 return 216;
2219 case MSH_HEX_343:
2220 if(name) *name = "Hexahedron 343";
2221 return 343;
2222 case MSH_HEX_512:
2223 if(name) *name = "Hexahedron 512";
2224 return 512;
2225 case MSH_HEX_729:
2226 if(name) *name = "Hexahedron 729";
2227 return 729;
2228 case MSH_HEX_1000:
2229 if(name) *name = "Hexahedron 1000";
2230 return 1000;
2231 case MSH_HEX_32:
2232 if(name) *name = "Hexahedron 32";
2233 return 8 + 12 * 2;
2234 case MSH_HEX_44:
2235 if(name) *name = "Hexahedron 44";
2236 return 8 + 12 * 3;
2237 case MSH_HEX_56:
2238 if(name) *name = "Hexahedron 56";
2239 return 8 + 12 * 4;
2240 case MSH_HEX_68:
2241 if(name) *name = "Hexahedron 68";
2242 return 8 + 12 * 5;
2243 case MSH_HEX_80:
2244 if(name) *name = "Hexahedron 80";
2245 return 8 + 12 * 6;
2246 case MSH_HEX_92:
2247 if(name) *name = "Hexahedron 92";
2248 return 8 + 12 * 7;
2249 case MSH_HEX_104:
2250 if(name) *name = "Hexahedron 104";
2251 return 8 + 12 * 8;
2252 case MSH_PRI_1:
2253 if(name) *name = "Prism 1";
2254 return 1;
2255 case MSH_PRI_6:
2256 if(name) *name = "Prism 6";
2257 return 6;
2258 case MSH_PRI_15:
2259 if(name) *name = "Prism 15";
2260 return 6 + 9;
2261 case MSH_PRI_18:
2262 if(name) *name = "Prism 18";
2263 return 6 + 9 + 3;
2264 case MSH_PRI_40:
2265 if(name) *name = "Prism 40";
2266 return 6 + 18 + 12 + 2 + 2 * 1;
2267 case MSH_PRI_75:
2268 if(name) *name = "Prism 75";
2269 return 6 + 27 + 27 + 6 + 3 * 3;
2270 case MSH_PRI_126:
2271 if(name) *name = "Prism 126";
2272 return 6 + 36 + 48 + 12 + 4 * 6;
2273 case MSH_PRI_196:
2274 if(name) *name = "Prism 196";
2275 return 6 + 45 + 75 + 20 + 5 * 10;
2276 case MSH_PRI_288:
2277 if(name) *name = "Prism 288";
2278 return 6 + 54 + 108 + 30 + 6 * 15;
2279 case MSH_PRI_405:
2280 if(name) *name = "Prism 405";
2281 return 6 + 63 + 147 + 42 + 7 * 21;
2282 case MSH_PRI_550:
2283 if(name) *name = "Prism 550";
2284 return 6 + 72 + 192 + 56 + 8 * 28;
2285 case MSH_PRI_24:
2286 if(name) *name = "Prism 24";
2287 return 6 + 9 * 2;
2288 case MSH_PRI_33:
2289 if(name) *name = "Prism 33";
2290 return 6 + 9 * 3;
2291 case MSH_PRI_42:
2292 if(name) *name = "Prism 42";
2293 return 6 + 9 * 4;
2294 case MSH_PRI_51:
2295 if(name) *name = "Prism 51";
2296 return 6 + 9 * 5;
2297 case MSH_PRI_60:
2298 if(name) *name = "Prism 60";
2299 return 6 + 9 * 6;
2300 case MSH_PRI_69:
2301 if(name) *name = "Prism 69";
2302 return 6 + 9 * 7;
2303 case MSH_PRI_78:
2304 if(name) *name = "Prism 78";
2305 return 6 + 9 * 8;
2306 case MSH_PYR_1:
2307 if(name) *name = "Pyramid 1";
2308 return 1;
2309 case MSH_PYR_5:
2310 if(name) *name = "Pyramid 5";
2311 return 5;
2312 case MSH_PYR_13:
2313 if(name) *name = "Pyramid 13";
2314 return 5 + 8;
2315 case MSH_PYR_14:
2316 if(name) *name = "Pyramid 14";
2317 return 5 + 8 + 1;
2318 case MSH_PYR_30:
2319 if(name) *name = "Pyramid 30";
2320 return 5 + 8 * 2 + 4 * 1 + 1 * 4 + 1;
2321 case MSH_PYR_55:
2322 if(name) *name = "Pyramid 55";
2323 return 5 + 8 * 3 + 4 * 3 + 1 * 9 + 5;
2324 case MSH_PYR_91:
2325 if(name) *name = "Pyramid 91";
2326 return 5 + 8 * 4 + 4 * 6 + 1 * 16 + 14;
2327 case MSH_PYR_140:
2328 if(name) *name = "Pyramid 140";
2329 return 5 + 8 * 5 + 4 * 10 + 1 * 25 + 30;
2330 case MSH_PYR_204:
2331 if(name) *name = "Pyramid 204";
2332 return 5 + 8 * 6 + 4 * 15 + 1 * 36 + 55;
2333 case MSH_PYR_285:
2334 if(name) *name = "Pyramid 285";
2335 return 5 + 8 * 7 + 4 * 21 + 1 * 49 + 91;
2336 case MSH_PYR_385:
2337 if(name) *name = "Pyramid 385";
2338 return 5 + 8 * 8 + 4 * 28 + 1 * 64 + 140;
2339 case MSH_PYR_21:
2340 if(name) *name = "Pyramid 21";
2341 return 5 + 8 * 2;
2342 case MSH_PYR_29:
2343 if(name) *name = "Pyramid 29";
2344 return 5 + 8 * 3;
2345 case MSH_PYR_37:
2346 if(name) *name = "Pyramid 37";
2347 return 5 + 8 * 4;
2348 case MSH_PYR_45:
2349 if(name) *name = "Pyramid 45";
2350 return 5 + 8 * 5;
2351 case MSH_PYR_53:
2352 if(name) *name = "Pyramid 53";
2353 return 5 + 8 * 6;
2354 case MSH_PYR_61:
2355 if(name) *name = "Pyramid 61";
2356 return 5 + 8 * 7;
2357 case MSH_PYR_69:
2358 if(name) *name = "Pyramid 69";
2359 return 5 + 8 * 8;
2360 case MSH_TRIH_4:
2361 if(name) *name = "Trihedron 4";
2362 return 4;
2363 case MSH_POLYH_:
2364 if(name) *name = "Polyhedron";
2365 return 0;
2366 case MSH_PNT_SUB:
2367 if(name) *name = "Point Xfem";
2368 return 1;
2369 case MSH_LIN_SUB:
2370 if(name) *name = "Line Xfem";
2371 return 2;
2372 case MSH_TRI_SUB:
2373 if(name) *name = "Triangle Xfem";
2374 return 3;
2375 case MSH_TET_SUB:
2376 if(name) *name = "Tetrahedron Xfem";
2377 return 4;
2378 default:
2379 Msg::Error("Unknown type of element %d", typeMSH);
2380 if(name) *name = "Unknown";
2381 return -1;
2382 }
2383 }
2384
getName()2385 std::string MElement::getName()
2386 {
2387 const char *name;
2388 MElement::getInfoMSH(getTypeForMSH(), &name);
2389 return name;
2390 }
2391
getVerticesIdForMSH(std::vector<int> & verts)2392 void MElement::getVerticesIdForMSH(std::vector<int> &verts)
2393 {
2394 int n = getNumVerticesForMSH();
2395 verts.resize(n);
2396 for(int i = 0; i < n; i++) verts[i] = (int)getVertex(i)->getIndex();
2397 }
2398
copy(std::map<int,MVertex * > & vertexMap,std::map<MElement *,MElement * > & newParents,std::map<MElement *,MElement * > & newDomains)2399 MElement *MElement::copy(std::map<int, MVertex *> &vertexMap,
2400 std::map<MElement *, MElement *> &newParents,
2401 std::map<MElement *, MElement *> &newDomains)
2402 {
2403 if(newDomains.count(this)) return newDomains.find(this)->second;
2404 std::vector<MVertex *> vmv;
2405 int eType = getTypeForMSH();
2406 MElement *eParent = getParent();
2407 if(getNumChildren() == 0) {
2408 for(std::size_t i = 0; i < getNumVertices(); i++) {
2409 MVertex *v = getVertex(i);
2410 std::size_t numV = v->getNum();
2411 if(vertexMap.count(numV))
2412 vmv.push_back(vertexMap[numV]);
2413 else {
2414 MVertex *mv = new MVertex(v->x(), v->y(), v->z(), nullptr, numV);
2415 vmv.push_back(mv);
2416 vertexMap[numV] = mv;
2417 }
2418 }
2419 }
2420 else {
2421 for(int i = 0; i < getNumChildren(); i++) {
2422 for(std::size_t j = 0; j < getChild(i)->getNumVertices(); j++) {
2423 MVertex *v = getChild(i)->getVertex(j);
2424 std::size_t numV = v->getNum();
2425 if(vertexMap.count(numV))
2426 vmv.push_back(vertexMap[numV]);
2427 else {
2428 MVertex *mv = new MVertex(v->x(), v->y(), v->z(), nullptr, numV);
2429 vmv.push_back(mv);
2430 vertexMap[numV] = mv;
2431 }
2432 }
2433 }
2434 }
2435
2436 MElement *parent = nullptr;
2437 if(eParent && !getDomain(0) && !getDomain(1)) {
2438 auto it = newParents.find(eParent);
2439 MElement *newParent;
2440 if(it == newParents.end()) {
2441 newParent = eParent->copy(vertexMap, newParents, newDomains);
2442 newParents[eParent] = newParent;
2443 }
2444 else
2445 newParent = it->second;
2446 parent = newParent;
2447 }
2448
2449 MElementFactory f;
2450 MElement *newEl =
2451 f.create(eType, vmv, getNum(), _partition, ownsParent(), 0, parent);
2452
2453 for(int i = 0; i < 2; i++) {
2454 MElement *dom = getDomain(i);
2455 if(!dom) continue;
2456 auto it = newDomains.find(dom);
2457 MElement *newDom;
2458 if(it == newDomains.end()) {
2459 newDom = dom->copy(vertexMap, newParents, newDomains);
2460 newDomains[dom] = newDom;
2461 }
2462 else
2463 newDom = newDomains.find(dom)->second;
2464 newEl->setDomain(newDom, i);
2465 }
2466 return newEl;
2467 }
2468
create(int type,std::vector<MVertex * > & v,std::size_t num,int part,bool owner,int parent,MElement * parent_ptr,MElement * d1,MElement * d2)2469 MElement *MElementFactory::create(int type, std::vector<MVertex *> &v,
2470 std::size_t num, int part, bool owner,
2471 int parent, MElement *parent_ptr,
2472 MElement *d1, MElement *d2)
2473 {
2474 switch(type) {
2475 case MSH_PNT: return new MPoint(v, num, part);
2476 case MSH_LIN_2: return new MLine(v, num, part);
2477 case MSH_LIN_3: return new MLine3(v, num, part);
2478 case MSH_LIN_4: return new MLineN(v, num, part);
2479 case MSH_LIN_5: return new MLineN(v, num, part);
2480 case MSH_LIN_6: return new MLineN(v, num, part);
2481 case MSH_LIN_7: return new MLineN(v, num, part);
2482 case MSH_LIN_8: return new MLineN(v, num, part);
2483 case MSH_LIN_9: return new MLineN(v, num, part);
2484 case MSH_LIN_10: return new MLineN(v, num, part);
2485 case MSH_LIN_11: return new MLineN(v, num, part);
2486 case MSH_LIN_B: return new MLineBorder(v, num, part, d1, d2);
2487 case MSH_LIN_C: return new MLineChild(v, num, part, owner, parent_ptr);
2488 case MSH_TRI_3: return new MTriangle(v, num, part);
2489 case MSH_TRI_6: return new MTriangle6(v, num, part);
2490 case MSH_TRI_10: return new MTriangleN(v, 3, num, part);
2491 case MSH_TRI_15: return new MTriangleN(v, 4, num, part);
2492 case MSH_TRI_21: return new MTriangleN(v, 5, num, part);
2493 case MSH_TRI_28: return new MTriangleN(v, 6, num, part);
2494 case MSH_TRI_36: return new MTriangleN(v, 7, num, part);
2495 case MSH_TRI_45: return new MTriangleN(v, 8, num, part);
2496 case MSH_TRI_55: return new MTriangleN(v, 9, num, part);
2497 case MSH_TRI_66: return new MTriangleN(v, 10, num, part);
2498 case MSH_TRI_9: return new MTriangleN(v, 3, num, part);
2499 case MSH_TRI_12: return new MTriangleN(v, 4, num, part);
2500 case MSH_TRI_15I: return new MTriangleN(v, 5, num, part);
2501 case MSH_TRI_18: return new MTriangleN(v, 6, num, part);
2502 case MSH_TRI_21I: return new MTriangleN(v, 7, num, part);
2503 case MSH_TRI_24: return new MTriangleN(v, 8, num, part);
2504 case MSH_TRI_27: return new MTriangleN(v, 9, num, part);
2505 case MSH_TRI_30: return new MTriangleN(v, 10, num, part);
2506 case MSH_TRI_B: return new MTriangleBorder(v, num, part, d1, d2);
2507 case MSH_QUA_4: return new MQuadrangle(v, num, part);
2508 case MSH_QUA_9: return new MQuadrangle9(v, num, part);
2509 case MSH_QUA_16: return new MQuadrangleN(v, 3, num, part);
2510 case MSH_QUA_25: return new MQuadrangleN(v, 4, num, part);
2511 case MSH_QUA_36: return new MQuadrangleN(v, 5, num, part);
2512 case MSH_QUA_49: return new MQuadrangleN(v, 6, num, part);
2513 case MSH_QUA_64: return new MQuadrangleN(v, 7, num, part);
2514 case MSH_QUA_81: return new MQuadrangleN(v, 8, num, part);
2515 case MSH_QUA_100: return new MQuadrangleN(v, 9, num, part);
2516 case MSH_QUA_121: return new MQuadrangleN(v, 10, num, part);
2517 case MSH_QUA_8: return new MQuadrangle8(v, num, part);
2518 case MSH_QUA_12: return new MQuadrangleN(v, 3, num, part);
2519 case MSH_QUA_16I: return new MQuadrangleN(v, 4, num, part);
2520 case MSH_QUA_20: return new MQuadrangleN(v, 5, num, part);
2521 case MSH_QUA_24: return new MQuadrangleN(v, 6, num, part);
2522 case MSH_QUA_28: return new MQuadrangleN(v, 7, num, part);
2523 case MSH_QUA_32: return new MQuadrangleN(v, 8, num, part);
2524 case MSH_QUA_36I: return new MQuadrangleN(v, 9, num, part);
2525 case MSH_QUA_40: return new MQuadrangleN(v, 10, num, part);
2526 case MSH_POLYG_: return new MPolygon(v, num, part, owner, parent_ptr);
2527 case MSH_POLYG_B: return new MPolygonBorder(v, num, part, d1, d2);
2528 case MSH_TET_4: return new MTetrahedron(v, num, part);
2529 case MSH_TET_10: return new MTetrahedron10(v, num, part);
2530 case MSH_HEX_8: return new MHexahedron(v, num, part);
2531 case MSH_HEX_20: return new MHexahedron20(v, num, part);
2532 case MSH_HEX_27: return new MHexahedron27(v, num, part);
2533 case MSH_PRI_6: return new MPrism(v, num, part);
2534 case MSH_PRI_15: return new MPrism15(v, num, part);
2535 case MSH_PRI_18: return new MPrism18(v, num, part);
2536 case MSH_PRI_40: return new MPrismN(v, 3, num, part);
2537 case MSH_PRI_75: return new MPrismN(v, 4, num, part);
2538 case MSH_PRI_126: return new MPrismN(v, 5, num, part);
2539 case MSH_PRI_196: return new MPrismN(v, 6, num, part);
2540 case MSH_PRI_288: return new MPrismN(v, 7, num, part);
2541 case MSH_PRI_405: return new MPrismN(v, 8, num, part);
2542 case MSH_PRI_550: return new MPrismN(v, 9, num, part);
2543 case MSH_PRI_24: return new MPrismN(v, 3, num, part);
2544 case MSH_PRI_33: return new MPrismN(v, 4, num, part);
2545 case MSH_PRI_42: return new MPrismN(v, 5, num, part);
2546 case MSH_PRI_51: return new MPrismN(v, 6, num, part);
2547 case MSH_PRI_60: return new MPrismN(v, 7, num, part);
2548 case MSH_PRI_69: return new MPrismN(v, 8, num, part);
2549 case MSH_PRI_78: return new MPrismN(v, 9, num, part);
2550 case MSH_PRI_1: return new MPrismN(v, 0, num, part);
2551 case MSH_TET_20: return new MTetrahedronN(v, 3, num, part);
2552 case MSH_TET_35: return new MTetrahedronN(v, 4, num, part);
2553 case MSH_TET_56: return new MTetrahedronN(v, 5, num, part);
2554 case MSH_TET_84: return new MTetrahedronN(v, 6, num, part);
2555 case MSH_TET_120: return new MTetrahedronN(v, 7, num, part);
2556 case MSH_TET_165: return new MTetrahedronN(v, 8, num, part);
2557 case MSH_TET_220: return new MTetrahedronN(v, 9, num, part);
2558 case MSH_TET_286: return new MTetrahedronN(v, 10, num, part);
2559 case MSH_TET_16: return new MTetrahedronN(v, 3, num, part);
2560 case MSH_TET_22: return new MTetrahedronN(v, 4, num, part);
2561 case MSH_TET_28: return new MTetrahedronN(v, 5, num, part);
2562 case MSH_TET_34: return new MTetrahedronN(v, 6, num, part);
2563 case MSH_TET_40: return new MTetrahedronN(v, 7, num, part);
2564 case MSH_TET_46: return new MTetrahedronN(v, 8, num, part);
2565 case MSH_TET_52: return new MTetrahedronN(v, 9, num, part);
2566 case MSH_TET_58: return new MTetrahedronN(v, 10, num, part);
2567 case MSH_POLYH_: return new MPolyhedron(v, num, part, owner, parent_ptr);
2568 case MSH_HEX_32: return new MHexahedronN(v, 3, num, part);
2569 case MSH_HEX_64: return new MHexahedronN(v, 3, num, part);
2570 case MSH_HEX_125: return new MHexahedronN(v, 4, num, part);
2571 case MSH_HEX_216: return new MHexahedronN(v, 5, num, part);
2572 case MSH_HEX_343: return new MHexahedronN(v, 6, num, part);
2573 case MSH_HEX_512: return new MHexahedronN(v, 7, num, part);
2574 case MSH_HEX_729: return new MHexahedronN(v, 8, num, part);
2575 case MSH_HEX_1000: return new MHexahedronN(v, 9, num, part);
2576 case MSH_PNT_SUB:
2577 return (parent_ptr) ? new MSubPoint(v, num, part, owner, parent_ptr) :
2578 new MSubPoint(v, num, part, owner, parent);
2579 case MSH_LIN_SUB:
2580 return (parent_ptr) ? new MSubLine(v, num, part, owner, parent_ptr) :
2581 new MSubLine(v, num, part, owner, parent);
2582 case MSH_TRI_SUB:
2583 return (parent_ptr) ? new MSubTriangle(v, num, part, owner, parent_ptr) :
2584 new MSubTriangle(v, num, part, owner, parent);
2585 case MSH_TET_SUB:
2586 return (parent_ptr) ? new MSubTetrahedron(v, num, part, owner, parent_ptr) :
2587 new MSubTetrahedron(v, num, part, owner, parent);
2588 case MSH_PYR_5: return new MPyramid(v, num, part);
2589 case MSH_PYR_13: return new MPyramidN(v, 2, num, part);
2590 case MSH_PYR_14: return new MPyramidN(v, 2, num, part);
2591 case MSH_PYR_30: return new MPyramidN(v, 3, num, part);
2592 case MSH_PYR_55: return new MPyramidN(v, 4, num, part);
2593 case MSH_PYR_91: return new MPyramidN(v, 5, num, part);
2594 case MSH_PYR_140: return new MPyramidN(v, 6, num, part);
2595 case MSH_PYR_204: return new MPyramidN(v, 7, num, part);
2596 case MSH_PYR_285: return new MPyramidN(v, 8, num, part);
2597 case MSH_PYR_385: return new MPyramidN(v, 9, num, part);
2598 case MSH_TRIH_4: return new MTrihedron(v, num, part);
2599 default: return nullptr;
2600 }
2601 }
2602
create(int num,int type,const std::vector<int> & data,GModel * model)2603 MElement *MElementFactory::create(int num, int type,
2604 const std::vector<int> &data, GModel *model)
2605 {
2606 // This should be rewritten: each element should register itself in a static
2607 // factory owned e.g. directly by MElement, and interpret its data by
2608 // itself. This would remove the ugly switch in the routine above.
2609
2610 int numVertices = MElement::getInfoMSH(type), startVertices = 0;
2611 if(data.size() && !numVertices) {
2612 startVertices = 1;
2613 numVertices = data[0];
2614 }
2615
2616 std::vector<MVertex *> vertices(numVertices);
2617 if((int)data.size() > startVertices + numVertices - 1) {
2618 for(int i = 0; i < numVertices; i++) {
2619 int numVertex = data[startVertices + i];
2620 MVertex *v = model->getMeshVertexByTag(numVertex);
2621 if(v) { vertices[i] = v; }
2622 else {
2623 Msg::Error("Unknown node %d in element %d", numVertex, num);
2624 return nullptr;
2625 }
2626 }
2627 }
2628 else {
2629 Msg::Error("Missing data in element %d", num);
2630 return nullptr;
2631 }
2632
2633 unsigned int part = 0;
2634 int startPartitions = startVertices + numVertices;
2635
2636 int parent = 0;
2637 if((type == MSH_PNT_SUB || type == MSH_LIN_SUB || type == MSH_TRI_SUB ||
2638 type == MSH_TET_SUB)) {
2639 parent = data[startPartitions];
2640 startPartitions += 1;
2641 }
2642
2643 std::vector<short> ghosts;
2644 if((int)data.size() > startPartitions) {
2645 int numPartitions = data[startPartitions];
2646 if(numPartitions > 0 &&
2647 (int)data.size() > startPartitions + numPartitions - 1) {
2648 part = data[startPartitions + 1];
2649 for(int i = 1; i < numPartitions; i++)
2650 ghosts.push_back(data[startPartitions + 1 + i]);
2651 }
2652 }
2653
2654 MElement *element = create(type, vertices, num, part, false, parent);
2655
2656 //for(std::size_t j = 0; j < ghosts.size(); j++) {
2657 // model->getGhostCells().insert(std::make_pair(element, ghosts[j]));
2658 //}
2659 if(part > model->getNumPartitions()) model->setNumPartitions(part);
2660
2661 return element;
2662 }
2663
skewness()2664 double MElement::skewness()
2665 {
2666 double minsk = 1.0;
2667 for(int i = 0; i < getNumFaces(); i++) {
2668 MFace f = getFace(i);
2669 if(f.getNumVertices() == 3) {
2670 // MTriangle t (f.getVertex(0),f.getVertex(1),f.getVertex(2));
2671 // minsk = std::min(minsk, t.etaShapeMeasure ());
2672 }
2673 else if(f.getNumVertices() == 4) {
2674 MQuadrangle q(f.getVertex(0), f.getVertex(1), f.getVertex(2),
2675 f.getVertex(3));
2676 minsk = std::min(minsk, q.etaShapeMeasure());
2677 }
2678 }
2679 return minsk;
2680 }
2681