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 "JacobianBasis.h"
7 #include "pointsGenerators.h"
8 #include "nodalBasis.h"
9 #include "BasisFactory.h"
10 #include "Numeric.h"
11 #include <cmath>
12
13 namespace {
14
15 template <class T>
calcMapFromIdealElement(int type,T & dSVec_dX,T & dSVec_dY,T & dSVec_dZ)16 void calcMapFromIdealElement(int type, T &dSVec_dX, T &dSVec_dY, T &dSVec_dZ)
17 {
18 // 2D scaling
19 switch(type) {
20 case TYPE_QUA:
21 case TYPE_HEX:
22 case TYPE_PYR: { // Quad, hex, pyramid -> square with side of length 1
23 dSVec_dX.scale(2.);
24 dSVec_dY.scale(2.);
25 break;
26 }
27 default: { // Tri, tet, prism: equilateral tri with side of length 1
28 static const double cTri[2] = {-1. / std::sqrt(3.), 2. / std::sqrt(3.)};
29 dSVec_dY.scale(cTri[1]);
30 dSVec_dY.axpy(dSVec_dX, cTri[0]);
31 break;
32 }
33 }
34
35 // 3D scaling
36 switch(type) {
37 case TYPE_HEX:
38 case TYPE_PRI: { // Hex, prism -> side of length 1 in z
39 dSVec_dZ.scale(2.);
40 break;
41 }
42 case TYPE_PYR: { // Pyramid -> height sqrt(2.)/2
43 static const double cPyr = sqrt(2.);
44 dSVec_dZ.scale(cPyr);
45 break;
46 }
47 case TYPE_TET: { // Tet: take into account (x, y) scaling to obtain regular
48 // tet
49 static const double cTet[3] = {-3. / 2 / std::sqrt(6.),
50 -1. / 2 / std::sqrt(2.), std::sqrt(1.5)};
51 dSVec_dZ.scale(cTet[2]);
52 dSVec_dZ.axpy(dSVec_dX, cTet[0]);
53 dSVec_dZ.axpy(dSVec_dY, cTet[1]);
54 break;
55 }
56 }
57 }
58
59 // Compute the determinant of a 3x3 matrix
calcDet3D(double M11,double M12,double M13,double M21,double M22,double M23,double M31,double M32,double M33)60 inline double calcDet3D(double M11, double M12, double M13, double M21,
61 double M22, double M23, double M31, double M32,
62 double M33)
63 {
64 return M11 * (M22 * M33 - M23 * M32) - M12 * (M21 * M33 - M23 * M31) +
65 M13 * (M21 * M32 - M22 * M31);
66 }
67
68 // Compute (signed) Jacobian determinant, and its partial derivative w.r.t.
69 // nodes coordinates, at sampling point 'i' of a 1D element (for which columns
70 // 2 and 3 of the Jacobian matrix are constant).
calcJDJ1D(double dxdX,double dxdY,double dxdZ,double dydX,double dydY,double dydZ,double dzdX,double dzdY,double dzdZ,int i,int numMapNodes,const fullMatrix<double> & dShapeMat_dX,fullMatrix<double> & JDJ)71 inline void calcJDJ1D(double dxdX, double dxdY, double dxdZ, double dydX,
72 double dydY, double dydZ, double dzdX, double dzdY,
73 double dzdZ, int i, int numMapNodes,
74 const fullMatrix<double> &dShapeMat_dX,
75 fullMatrix<double> &JDJ)
76 {
77 for(int j = 0; j < numMapNodes; j++) {
78 const double &dPhidX = dShapeMat_dX(i, j);
79 JDJ(i, j) = dPhidX * (dydY * dzdZ - dzdY * dydZ);
80 JDJ(i, j + numMapNodes) = dPhidX * (dzdY * dxdZ - dxdY * dzdZ);
81 JDJ(i, j + 2 * numMapNodes) = dPhidX * (dxdY * dydZ - dydY * dxdZ);
82 }
83 JDJ(i, 3 * numMapNodes) =
84 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
85 }
86
87 // Compute (signed) Jacobian determinant, and its partial derivative w.r.t.
88 // nodes coordinates, at sampling point 'i' of a 2D element (for which column
89 // 3 of the Jacobian matrix is constant).
calcJDJ2D(double dxdX,double dxdY,double dxdZ,double dydX,double dydY,double dydZ,double dzdX,double dzdY,double dzdZ,int i,int numMapNodes,const fullMatrix<double> & dShapeMat_dX,const fullMatrix<double> & dShapeMat_dY,fullMatrix<double> & JDJ)90 inline void calcJDJ2D(double dxdX, double dxdY, double dxdZ, double dydX,
91 double dydY, double dydZ, double dzdX, double dzdY,
92 double dzdZ, int i, int numMapNodes,
93 const fullMatrix<double> &dShapeMat_dX,
94 const fullMatrix<double> &dShapeMat_dY,
95 fullMatrix<double> &JDJ)
96 {
97 for(int j = 0; j < numMapNodes; j++) {
98 const double &dPhidX = dShapeMat_dX(i, j);
99 const double &dPhidY = dShapeMat_dY(i, j);
100 JDJ(i, j) = dPhidX * (dydY * dzdZ - dzdY * dydZ) +
101 dPhidY * (dzdX * dydZ - dydX * dzdZ);
102 JDJ(i, j + numMapNodes) = dPhidX * (dzdY * dxdZ - dxdY * dzdZ) +
103 dPhidY * (dxdX * dzdZ - dzdX * dxdZ);
104 JDJ(i, j + 2 * numMapNodes) = dPhidX * (dxdY * dydZ - dydY * dxdZ) +
105 dPhidY * (dydX * dxdZ - dxdX * dydZ);
106 }
107 JDJ(i, 3 * numMapNodes) =
108 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
109 }
110
111 // Compute (signed) Jacobian determinant, and its partial derivative w.r.t.
112 // nodes coordinates, at sampling point 'i' of a 3D element
calcJDJ3D(double dxdX,double dxdY,double dxdZ,double dydX,double dydY,double dydZ,double dzdX,double dzdY,double dzdZ,int i,int numMapNodes,const fullMatrix<double> & dShapeMat_dX,const fullMatrix<double> & dShapeMat_dY,const fullMatrix<double> & dShapeMat_dZ,fullMatrix<double> & JDJ)113 inline void calcJDJ3D(double dxdX, double dxdY, double dxdZ, double dydX,
114 double dydY, double dydZ, double dzdX, double dzdY,
115 double dzdZ, int i, int numMapNodes,
116 const fullMatrix<double> &dShapeMat_dX,
117 const fullMatrix<double> &dShapeMat_dY,
118 const fullMatrix<double> &dShapeMat_dZ,
119 fullMatrix<double> &JDJ)
120 {
121 for(int j = 0; j < numMapNodes; j++) {
122 const double &dPhidX = dShapeMat_dX(i, j);
123 const double &dPhidY = dShapeMat_dY(i, j);
124 const double &dPhidZ = dShapeMat_dZ(i, j);
125 JDJ(i, j) = dPhidX * (dydY * dzdZ - dzdY * dydZ) +
126 dPhidY * (dzdX * dydZ - dydX * dzdZ) +
127 dPhidZ * (dydX * dzdY - dzdX * dydY);
128 JDJ(i, j + numMapNodes) = dPhidX * (dzdY * dxdZ - dxdY * dzdZ) +
129 dPhidY * (dxdX * dzdZ - dzdX * dxdZ) +
130 dPhidZ * (dzdX * dxdY - dxdX * dzdY);
131 JDJ(i, j + 2 * numMapNodes) = dPhidX * (dxdY * dydZ - dydY * dxdZ) +
132 dPhidY * (dydX * dxdZ - dxdX * dydZ) +
133 dPhidZ * (dxdX * dydY - dydX * dxdY);
134 }
135 JDJ(i, 3 * numMapNodes) =
136 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
137 }
138
139 } // namespace
140
141 // GradientBasis contains the information for computing the partial
142 // derivatives of the mapping components w.r.t. to each reference coordinate
143 // for a given element type.
144 // Those partial derivatives are the components of the Jacobian matrix,
145 // and the gradients are the rows of the Jacobian matrix.
GradientBasis(int elementTag,FuncSpaceData fsdata)146 GradientBasis::GradientBasis(int elementTag, FuncSpaceData fsdata)
147 : _elementTag(elementTag), _data(fsdata)
148 {
149 // Matrix dShapeMat_dX, when multiplied by Lagrange coefficients
150 // (i.e. node coordinates), gives the derivative of the mapping
151 // w.r.t. the first reference coordinate at some sampling points.
152 // Those sampling points is determined by 'fsdata'.
153 // The ordering of the sampling points is "ordered" (see pointsGenerator.cpp)
154 // and is thus different from the Gmsh ordering convention. This order is
155 // useful for efficiently converting Lagrange to Bezier coefficients.
156 fullMatrix<double> samplingPoints;
157 gmshGenerateOrderedPoints(fsdata, samplingPoints);
158 const int numSampPnts = samplingPoints.size1();
159
160 // Store partial derivatives of shape functions at sampling points
161 fullMatrix<double> allDPsi;
162 const nodalBasis *mapBasis = BasisFactory::getNodalBasis(_elementTag);
163 mapBasis->df(samplingPoints, allDPsi);
164 const int numMapNodes = allDPsi.size2();
165
166 dShapeMat_dX.resize(numSampPnts, numMapNodes);
167 dShapeMat_dY.resize(numSampPnts, numMapNodes);
168 dShapeMat_dZ.resize(numSampPnts, numMapNodes);
169 for(int i = 0; i < numSampPnts; i++) {
170 for(int j = 0; j < numMapNodes; j++) {
171 dShapeMat_dX(i, j) = allDPsi(3 * i + 0, j);
172 dShapeMat_dY(i, j) = allDPsi(3 * i + 1, j);
173 dShapeMat_dZ(i, j) = allDPsi(3 * i + 2, j);
174 }
175 }
176
177 dShapeIdealMat_dX = dShapeMat_dX;
178 dShapeIdealMat_dY = dShapeMat_dY;
179 dShapeIdealMat_dZ = dShapeMat_dZ;
180 mapFromIdealElement(_data.getType(), dShapeIdealMat_dX, dShapeIdealMat_dY,
181 dShapeIdealMat_dZ);
182 }
183
getIdealGradientsFromNodes(const fullMatrix<double> & nodesCoord,fullMatrix<double> * dxyzdX,fullMatrix<double> * dxyzdY,fullMatrix<double> * dxyzdZ) const184 void GradientBasis::getIdealGradientsFromNodes(
185 const fullMatrix<double> &nodesCoord, fullMatrix<double> *dxyzdX,
186 fullMatrix<double> *dxyzdY, fullMatrix<double> *dxyzdZ) const
187 {
188 if(dxyzdX) dShapeIdealMat_dX.mult(nodesCoord, *dxyzdX);
189 if(dxyzdY) dShapeIdealMat_dY.mult(nodesCoord, *dxyzdY);
190 if(dxyzdZ) dShapeIdealMat_dZ.mult(nodesCoord, *dxyzdZ);
191 }
192
getGradientsFromNodes(const fullMatrix<double> & nodesCoord,fullMatrix<double> * dxyzdX,fullMatrix<double> * dxyzdY,fullMatrix<double> * dxyzdZ) const193 void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodesCoord,
194 fullMatrix<double> *dxyzdX,
195 fullMatrix<double> *dxyzdY,
196 fullMatrix<double> *dxyzdZ) const
197 {
198 if(dxyzdX) dShapeMat_dX.mult(nodesCoord, *dxyzdX);
199 if(dxyzdY) dShapeMat_dY.mult(nodesCoord, *dxyzdY);
200 if(dxyzdZ) dShapeMat_dZ.mult(nodesCoord, *dxyzdZ);
201 }
202
getAllGradientsFromNodes(const fullMatrix<double> & nodesCoord,fullMatrix<double> & dxyzdXYZ) const203 void GradientBasis::getAllGradientsFromNodes(
204 const fullMatrix<double> &nodesCoord, fullMatrix<double> &dxyzdXYZ) const
205 {
206 fullMatrix<double> prox;
207 prox.setAsProxy(dxyzdXYZ, 0, 3);
208 dShapeMat_dX.mult(nodesCoord, prox);
209
210 prox.setAsProxy(dxyzdXYZ, 3, 3);
211 dShapeMat_dY.mult(nodesCoord, prox);
212
213 prox.setAsProxy(dxyzdXYZ, 6, 3);
214 dShapeMat_dZ.mult(nodesCoord, prox);
215 }
216
getAllIdealGradientsFromNodes(const fullMatrix<double> & nodesCoord,fullMatrix<double> & dxyzdXYZ) const217 void GradientBasis::getAllIdealGradientsFromNodes(
218 const fullMatrix<double> &nodesCoord, fullMatrix<double> &dxyzdXYZ) const
219 {
220 fullMatrix<double> prox;
221 prox.setAsProxy(dxyzdXYZ, 0, 3);
222 dShapeIdealMat_dX.mult(nodesCoord, prox);
223
224 prox.setAsProxy(dxyzdXYZ, 3, 3);
225 dShapeIdealMat_dY.mult(nodesCoord, prox);
226
227 prox.setAsProxy(dxyzdXYZ, 6, 3);
228 dShapeIdealMat_dZ.mult(nodesCoord, prox);
229 }
230
mapFromIdealElement(int type,fullMatrix<double> & dSMat_dX,fullMatrix<double> & dSMat_dY,fullMatrix<double> & dSMat_dZ)231 void GradientBasis::mapFromIdealElement(int type, fullMatrix<double> &dSMat_dX,
232 fullMatrix<double> &dSMat_dY,
233 fullMatrix<double> &dSMat_dZ)
234 {
235 calcMapFromIdealElement(type, dSMat_dX, dSMat_dY, dSMat_dZ);
236 }
237
mapFromIdealElement(int type,fullVector<double> & dSVec_dX,fullVector<double> & dSVec_dY,fullVector<double> & dSVec_dZ)238 void GradientBasis::mapFromIdealElement(int type, fullVector<double> &dSVec_dX,
239 fullVector<double> &dSVec_dY,
240 fullVector<double> &dSVec_dZ)
241 {
242 calcMapFromIdealElement(type, dSVec_dX, dSVec_dY, dSVec_dZ);
243 }
244
mapFromIdealElement(int type,double jac[3][3])245 void GradientBasis::mapFromIdealElement(int type, double jac[3][3])
246 {
247 fullMatrix<double> dxyzdX(jac[0], 1, 3), dxyzdY(jac[1], 1, 3),
248 dxyzdZ(jac[2], 1, 3);
249 mapFromIdealElement(type, dxyzdX, dxyzdY, dxyzdZ);
250 }
251
JacobianBasis(int elementTag,FuncSpaceData data)252 JacobianBasis::JacobianBasis(int elementTag, FuncSpaceData data)
253 : _elementTag(elementTag), _data(data), _dim(data.getDimension())
254 {
255 const int parentType = data.getType();
256 const int primJacobianOrder = jacobianOrder(parentType, 1);
257
258 fullMatrix<double> samplingPoints;
259 gmshGeneratePoints(data, samplingPoints);
260 numSamplingPnts = samplingPoints.size1();
261
262 // Store shape function gradients of mapping at Jacobian nodes
263 _gradBasis = BasisFactory::getGradientBasis(elementTag, data);
264
265 // Compute matrix for lifting from primary Jacobian basis to Jacobian basis
266 int primJacType = ElementType::getType(parentType, primJacobianOrder, false);
267 const nodalBasis *primJacBasis = BasisFactory::getNodalBasis(primJacType);
268 numPrimSamplingPnts = primJacBasis->getNumShapeFunctions();
269
270 matrixPrimJac2Jac.resize(numSamplingPnts, numPrimSamplingPnts);
271 primJacBasis->f(samplingPoints, matrixPrimJac2Jac);
272
273 // Compute shape function gradients of primary mapping at barycenter, in order
274 // to compute normal to straight element
275 const int primMapType = ElementType::getType(parentType, 1, false);
276 const nodalBasis *primMapBasis = BasisFactory::getNodalBasis(primMapType);
277 numPrimMapNodes = primMapBasis->getNumShapeFunctions();
278
279 double xBar = 0., yBar = 0., zBar = 0.;
280 double barycenter[3] = {0., 0., 0.};
281 for(int i = 0; i < numPrimMapNodes; i++) {
282 for(int j = 0; j < primMapBasis->points.size2(); ++j) {
283 barycenter[j] += primMapBasis->points(i, j);
284 }
285 }
286 barycenter[0] /= numPrimMapNodes;
287 barycenter[1] /= numPrimMapNodes;
288 barycenter[2] /= numPrimMapNodes;
289
290 double(*barDPsi)[3] = new double[numPrimMapNodes][3];
291 primMapBasis->df(xBar, yBar, zBar, barDPsi);
292
293 dPrimBaryShape_dX.resize(numPrimMapNodes);
294 dPrimBaryShape_dY.resize(numPrimMapNodes);
295 dPrimBaryShape_dZ.resize(numPrimMapNodes);
296 for(int j = 0; j < numPrimMapNodes; j++) {
297 dPrimBaryShape_dX(j) = barDPsi[j][0];
298 dPrimBaryShape_dY(j) = barDPsi[j][1];
299 dPrimBaryShape_dZ(j) = barDPsi[j][2];
300 }
301
302 dPrimBaryIdealShape_dX = dPrimBaryShape_dX;
303 dPrimBaryIdealShape_dY = dPrimBaryShape_dY;
304 dPrimBaryIdealShape_dZ = dPrimBaryShape_dZ;
305 _gradBasis->mapFromIdealElement(
306 dPrimBaryIdealShape_dX, dPrimBaryIdealShape_dY, dPrimBaryIdealShape_dZ);
307
308 delete[] barDPsi;
309
310 // Compute "fast" Jacobian evaluation matrices (at 1st order nodes +
311 // barycenter)
312 numSamplingPntsFast = numPrimMapNodes + 1;
313 fullMatrix<double> lagPointsFast(numSamplingPntsFast, 3); // Sampling points
314 lagPointsFast.copy(primMapBasis->points, 0, numPrimMapNodes, 0,
315 primMapBasis->points.size2(), 0, 0); // 1st order nodes
316 lagPointsFast(numPrimMapNodes, 0) = barycenter[0]; // Last point = barycenter
317 lagPointsFast(numPrimMapNodes, 1) = barycenter[1];
318 lagPointsFast(numPrimMapNodes, 2) = barycenter[2];
319
320 fullMatrix<double> allDPsiFast;
321 const nodalBasis *mapBasis = BasisFactory::getNodalBasis(_elementTag);
322 mapBasis->df(lagPointsFast, allDPsiFast);
323 numMapNodes = mapBasis->getNumShapeFunctions();
324
325 dFastShapeMat_dX.resize(numSamplingPntsFast, numMapNodes);
326 dFastShapeMat_dY.resize(numSamplingPntsFast, numMapNodes);
327 dFastShapeMat_dZ.resize(numSamplingPntsFast, numMapNodes);
328 for(int i = 0; i < numSamplingPntsFast; i++) {
329 for(int j = 0; j < numMapNodes; j++) {
330 dFastShapeMat_dX(i, j) = allDPsiFast(3 * i + 0, j);
331 dFastShapeMat_dY(i, j) = allDPsiFast(3 * i + 1, j);
332 dFastShapeMat_dZ(i, j) = allDPsiFast(3 * i + 2, j);
333 }
334 }
335 }
336
337 // Computes (unit) normals to straight line element at barycenter (with norm of
338 // gradient as return value)
getPrimNormals1D(const fullMatrix<double> & nodesXYZ,fullMatrix<double> & result) const339 double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ,
340 fullMatrix<double> &result) const
341 {
342 fullVector<double> dxyzdXbar(3);
343 for(int j = 0; j < numPrimMapNodes; j++) {
344 dxyzdXbar(0) += dPrimBaryShape_dX(j) * nodesXYZ(j, 0);
345 dxyzdXbar(1) += dPrimBaryShape_dX(j) * nodesXYZ(j, 1);
346 dxyzdXbar(2) += dPrimBaryShape_dX(j) * nodesXYZ(j, 2);
347 }
348
349 if((fabs(dxyzdXbar(0)) >= fabs(dxyzdXbar(1)) &&
350 fabs(dxyzdXbar(0)) >= fabs(dxyzdXbar(2))) ||
351 (fabs(dxyzdXbar(1)) >= fabs(dxyzdXbar(0)) &&
352 fabs(dxyzdXbar(1)) >= fabs(dxyzdXbar(2)))) {
353 result(0, 0) = dxyzdXbar(1);
354 result(0, 1) = -dxyzdXbar(0);
355 result(0, 2) = 0.;
356 }
357 else {
358 result(0, 0) = 0.;
359 result(0, 1) = dxyzdXbar(2);
360 result(0, 2) = -dxyzdXbar(1);
361 }
362 const double norm0 =
363 sqrt(result(0, 0) * result(0, 0) + result(0, 1) * result(0, 1) +
364 result(0, 2) * result(0, 2));
365 result(0, 0) /= norm0;
366 result(0, 1) /= norm0;
367 result(0, 2) /= norm0;
368
369 result(1, 2) = dxyzdXbar(0) * result(0, 1) - dxyzdXbar(1) * result(0, 0);
370 result(1, 1) = -dxyzdXbar(0) * result(0, 2) + dxyzdXbar(2) * result(0, 0);
371 result(1, 0) = dxyzdXbar(1) * result(0, 2) - dxyzdXbar(2) * result(0, 1);
372 const double norm1 =
373 sqrt(result(1, 0) * result(1, 0) + result(1, 1) * result(1, 1) +
374 result(1, 2) * result(1, 2));
375 result(1, 0) /= norm1;
376 result(1, 1) /= norm1;
377 result(1, 2) /= norm1;
378
379 return sqrt(dxyzdXbar(0) * dxyzdXbar(0) + dxyzdXbar(1) * dxyzdXbar(1) +
380 dxyzdXbar(2) * dxyzdXbar(2));
381 }
382
383 // Computes (unit) normal to straight surface element at barycenter (with norm
384 // as return value)
getPrimNormal2D(const fullMatrix<double> & nodesXYZ,fullMatrix<double> & result,bool ideal) const385 double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ,
386 fullMatrix<double> &result,
387 bool ideal) const
388 {
389 const fullVector<double> &gSX =
390 ideal ? dPrimBaryIdealShape_dX : dPrimBaryShape_dX;
391 const fullVector<double> &gSY =
392 ideal ? dPrimBaryIdealShape_dY : dPrimBaryShape_dY;
393 fullMatrix<double> dxyzdX(1, 3), dxyzdY(1, 3);
394 for(int j = 0; j < numPrimMapNodes; j++) {
395 dxyzdX(0, 0) += gSX(j) * nodesXYZ(j, 0);
396 dxyzdX(0, 1) += gSX(j) * nodesXYZ(j, 1);
397 dxyzdX(0, 2) += gSX(j) * nodesXYZ(j, 2);
398 dxyzdY(0, 0) += gSY(j) * nodesXYZ(j, 0);
399 dxyzdY(0, 1) += gSY(j) * nodesXYZ(j, 1);
400 dxyzdY(0, 2) += gSY(j) * nodesXYZ(j, 2);
401 }
402
403 result(0, 2) = dxyzdX(0, 0) * dxyzdY(0, 1) - dxyzdX(0, 1) * dxyzdY(0, 0);
404 result(0, 1) = -dxyzdX(0, 0) * dxyzdY(0, 2) + dxyzdX(0, 2) * dxyzdY(0, 0);
405 result(0, 0) = dxyzdX(0, 1) * dxyzdY(0, 2) - dxyzdX(0, 2) * dxyzdY(0, 1);
406 double norm0 =
407 sqrt(result(0, 0) * result(0, 0) + result(0, 1) * result(0, 1) +
408 result(0, 2) * result(0, 2));
409 const double invNorm0 = 1. / norm0;
410 result(0, 0) *= invNorm0;
411 result(0, 1) *= invNorm0;
412 result(0, 2) *= invNorm0;
413
414 return norm0;
415 }
416
417 // Returns absolute value of Jacobian of straight volume element at barycenter
getPrimJac3D(const fullMatrix<double> & nodesXYZ,bool ideal) const418 double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ,
419 bool ideal) const
420 {
421 const fullVector<double> &gSX =
422 ideal ? dPrimBaryIdealShape_dX : dPrimBaryShape_dX;
423 const fullVector<double> &gSY =
424 ideal ? dPrimBaryIdealShape_dY : dPrimBaryShape_dY;
425 const fullVector<double> &gSZ =
426 ideal ? dPrimBaryIdealShape_dZ : dPrimBaryShape_dZ;
427 fullMatrix<double> dxyzdX(1, 3), dxyzdY(1, 3), dxyzdZ(1, 3);
428 for(int j = 0; j < numPrimMapNodes; j++) {
429 dxyzdX(0, 0) += gSX(j) * nodesXYZ(j, 0);
430 dxyzdX(0, 1) += gSX(j) * nodesXYZ(j, 1);
431 dxyzdX(0, 2) += gSX(j) * nodesXYZ(j, 2);
432 dxyzdY(0, 0) += gSY(j) * nodesXYZ(j, 0);
433 dxyzdY(0, 1) += gSY(j) * nodesXYZ(j, 1);
434 dxyzdY(0, 2) += gSY(j) * nodesXYZ(j, 2);
435 dxyzdZ(0, 0) += gSZ(j) * nodesXYZ(j, 0);
436 dxyzdZ(0, 1) += gSZ(j) * nodesXYZ(j, 1);
437 dxyzdZ(0, 2) += gSZ(j) * nodesXYZ(j, 2);
438 }
439 double dJ = fabs(calcDet3D(dxyzdX(0, 0), dxyzdY(0, 0), dxyzdZ(0, 0),
440 dxyzdX(0, 1), dxyzdY(0, 1), dxyzdZ(0, 1),
441 dxyzdX(0, 2), dxyzdY(0, 2), dxyzdZ(0, 2)));
442 return dJ;
443 }
444
445 // Calculate (signed, possibly scaled) Jacobian for one element, with normal
446 // vectors to straight element for regularization. Evaluation points depend on
447 // the given matrices for shape function gradients.
getJacobianGeneral(int nSamplingPnts,const fullMatrix<double> & dSMat_dX,const fullMatrix<double> & dSMat_dY,const fullMatrix<double> & dSMat_dZ,const fullMatrix<double> & nodesXYZ,bool idealNorm,bool scaling,fullVector<double> & jacobian,const fullMatrix<double> * geomNormals) const448 void JacobianBasis::getJacobianGeneral(
449 int nSamplingPnts, const fullMatrix<double> &dSMat_dX,
450 const fullMatrix<double> &dSMat_dY, const fullMatrix<double> &dSMat_dZ,
451 const fullMatrix<double> &nodesXYZ, bool idealNorm, bool scaling,
452 fullVector<double> &jacobian, const fullMatrix<double> *geomNormals) const
453 {
454 switch(_dim) {
455 case 0: {
456 for(int i = 0; i < nSamplingPnts; i++) jacobian(i) = 1.;
457 } break;
458 case 1: {
459 fullMatrix<double> normals(2, 3);
460 const double invScale = getPrimNormals1D(nodesXYZ, normals);
461 if(geomNormals) normals.setAll(*geomNormals);
462 if(scaling) {
463 if(invScale == 0) {
464 for(int i = 0; i < nSamplingPnts; i++) jacobian(i) = 0;
465 return;
466 }
467 const double scale = 1. / invScale;
468 // Faster to scale 1 normal than afterwards
469 normals(0, 0) *= scale;
470 normals(0, 1) *= scale;
471 normals(0, 2) *= scale;
472 }
473 fullMatrix<double> dxyzdX(nSamplingPnts, 3);
474 dSMat_dX.mult(nodesXYZ, dxyzdX);
475 for(int i = 0; i < nSamplingPnts; i++) {
476 const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
477 &dzdX = dxyzdX(i, 2);
478 const double &dxdY = normals(0, 0), &dydY = normals(0, 1),
479 &dzdY = normals(0, 2);
480 const double &dxdZ = normals(1, 0), &dydZ = normals(1, 1),
481 &dzdZ = normals(1, 2);
482 jacobian(i) =
483 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
484 }
485 } break;
486 case 2: {
487 fullMatrix<double> normal(1, 3);
488 const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
489 if(geomNormals) normal.setAll(*geomNormals);
490 if(scaling) {
491 if(invScale == 0) {
492 for(int i = 0; i < nSamplingPnts; i++) jacobian(i) = 0;
493 return;
494 }
495 const double scale = 1. / invScale;
496 // Faster to scale normal than afterwards
497 normal(0, 0) *= scale;
498 normal(0, 1) *= scale;
499 normal(0, 2) *= scale;
500 }
501 fullMatrix<double> dxyzdX(nSamplingPnts, 3), dxyzdY(nSamplingPnts, 3);
502 dSMat_dX.mult(nodesXYZ, dxyzdX);
503 dSMat_dY.mult(nodesXYZ, dxyzdY);
504 for(int i = 0; i < nSamplingPnts; i++) {
505 const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
506 &dzdX = dxyzdX(i, 2);
507 const double &dxdY = dxyzdY(i, 0), &dydY = dxyzdY(i, 1),
508 &dzdY = dxyzdY(i, 2);
509 const double &dxdZ = normal(0, 0), &dydZ = normal(0, 1),
510 &dzdZ = normal(0, 2);
511 jacobian(i) =
512 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
513 }
514 } break;
515 case 3: {
516 fullMatrix<double> dum;
517 fullMatrix<double> dxyzdX(nSamplingPnts, 3), dxyzdY(nSamplingPnts, 3),
518 dxyzdZ(nSamplingPnts, 3);
519 dSMat_dX.mult(nodesXYZ, dxyzdX);
520 dSMat_dY.mult(nodesXYZ, dxyzdY);
521 dSMat_dZ.mult(nodesXYZ, dxyzdZ);
522 for(int i = 0; i < nSamplingPnts; i++) {
523 const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
524 &dzdX = dxyzdX(i, 2);
525 const double &dxdY = dxyzdY(i, 0), &dydY = dxyzdY(i, 1),
526 &dzdY = dxyzdY(i, 2);
527 const double &dxdZ = dxyzdZ(i, 0), &dydZ = dxyzdZ(i, 1),
528 &dzdZ = dxyzdZ(i, 2);
529 jacobian(i) =
530 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
531 }
532 if(scaling) {
533 const double scale = 1. / getPrimJac3D(nodesXYZ);
534 jacobian.scale(scale);
535 }
536 } break;
537 }
538 }
539
540 // Calculate (signed, possibly scaled) Jacobian for several elements, with
541 // normal vectors to straight elements for regularization. Evaluation points
542 // depend on the given matrices for shape function gradients. TODO: Optimize
543 // and test 1D & 2D
getJacobianGeneral(int nSamplingPnts,const fullMatrix<double> & dSMat_dX,const fullMatrix<double> & dSMat_dY,const fullMatrix<double> & dSMat_dZ,const fullMatrix<double> & nodesX,const fullMatrix<double> & nodesY,const fullMatrix<double> & nodesZ,bool idealNorm,bool scaling,fullMatrix<double> & jacobian,const fullMatrix<double> * geomNormals) const544 void JacobianBasis::getJacobianGeneral(
545 int nSamplingPnts, const fullMatrix<double> &dSMat_dX,
546 const fullMatrix<double> &dSMat_dY, const fullMatrix<double> &dSMat_dZ,
547 const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
548 const fullMatrix<double> &nodesZ, bool idealNorm, bool scaling,
549 fullMatrix<double> &jacobian, const fullMatrix<double> *geomNormals) const
550 {
551 switch(_dim) {
552 case 0: {
553 const int numEl = nodesX.size2();
554 for(int iEl = 0; iEl < numEl; iEl++)
555 for(int i = 0; i < nSamplingPnts; i++) jacobian(i, iEl) = 1.;
556 } break;
557 case 1: {
558 const int numEl = nodesX.size2();
559 fullMatrix<double> dxdX(nSamplingPnts, numEl);
560 fullMatrix<double> dydX(nSamplingPnts, numEl);
561 fullMatrix<double> dzdX(nSamplingPnts, numEl);
562 dSMat_dX.mult(nodesX, dxdX);
563 dSMat_dX.mult(nodesY, dydX);
564 dSMat_dX.mult(nodesZ, dzdX);
565 for(int iEl = 0; iEl < numEl; iEl++) {
566 fullMatrix<double> nodesXYZ(numPrimMapNodes, 3);
567 for(int i = 0; i < numPrimMapNodes; i++) {
568 nodesXYZ(i, 0) = nodesX(i, iEl);
569 nodesXYZ(i, 1) = nodesY(i, iEl);
570 nodesXYZ(i, 2) = nodesZ(i, iEl);
571 }
572 fullMatrix<double> normals(2, 3);
573 const double invScale = getPrimNormals1D(nodesXYZ, normals);
574 if(geomNormals) normals.setAll(*geomNormals);
575 if(scaling) {
576 if(invScale == 0) {
577 for(int i = 0; i < nSamplingPnts; i++) jacobian(i, iEl) = 0;
578 continue;
579 }
580 const double scale = 1. / invScale;
581 // Faster to scale 1 normal than afterwards
582 normals(0, 0) *= scale;
583 normals(0, 1) *= scale;
584 normals(0, 2) *= scale;
585 }
586 const double &dxdY = normals(0, 0), &dydY = normals(0, 1),
587 &dzdY = normals(0, 2);
588 const double &dxdZ = normals(1, 0), &dydZ = normals(1, 1),
589 &dzdZ = normals(1, 2);
590 for(int i = 0; i < nSamplingPnts; i++)
591 jacobian(i, iEl) = calcDet3D(dxdX(i, iEl), dxdY, dxdZ, dydX(i, iEl),
592 dydY, dydZ, dzdX(i, iEl), dzdY, dzdZ);
593 }
594 } break;
595 case 2: {
596 const int numEl = nodesX.size2();
597 fullMatrix<double> dxdX(nSamplingPnts, numEl);
598 fullMatrix<double> dydX(nSamplingPnts, numEl);
599 fullMatrix<double> dzdX(nSamplingPnts, numEl);
600 fullMatrix<double> dxdY(nSamplingPnts, numEl);
601 fullMatrix<double> dydY(nSamplingPnts, numEl);
602 fullMatrix<double> dzdY(nSamplingPnts, numEl);
603 dSMat_dX.mult(nodesX, dxdX);
604 dSMat_dX.mult(nodesY, dydX);
605 dSMat_dX.mult(nodesZ, dzdX);
606 dSMat_dY.mult(nodesX, dxdY);
607 dSMat_dY.mult(nodesY, dydY);
608 dSMat_dY.mult(nodesZ, dzdY);
609 for(int iEl = 0; iEl < numEl; iEl++) {
610 fullMatrix<double> nodesXYZ(numPrimMapNodes, 3);
611 for(int i = 0; i < numPrimMapNodes; i++) {
612 nodesXYZ(i, 0) = nodesX(i, iEl);
613 nodesXYZ(i, 1) = nodesY(i, iEl);
614 nodesXYZ(i, 2) = nodesZ(i, iEl);
615 }
616 fullMatrix<double> normal(1, 3);
617 const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
618 if(geomNormals) normal.setAll(*geomNormals);
619 if(scaling) {
620 if(invScale == 0) {
621 for(int i = 0; i < nSamplingPnts; i++) jacobian(i, iEl) = 0;
622 continue;
623 }
624 const double scale = 1. / invScale;
625 // Faster to scale normal than afterwards
626 normal(0, 0) *= scale;
627 normal(0, 1) *= scale;
628 normal(0, 2) *= scale;
629 }
630 const double &dxdZ = normal(0, 0), &dydZ = normal(0, 1),
631 &dzdZ = normal(0, 2);
632 for(int i = 0; i < nSamplingPnts; i++)
633 jacobian(i, iEl) =
634 calcDet3D(dxdX(i, iEl), dxdY(i, iEl), dxdZ, dydX(i, iEl),
635 dydY(i, iEl), dydZ, dzdX(i, iEl), dzdY(i, iEl), dzdZ);
636 }
637 } break;
638 case 3: {
639 const int numEl = nodesX.size2();
640 fullMatrix<double> dxdX(nSamplingPnts, numEl);
641 fullMatrix<double> dydX(nSamplingPnts, numEl);
642 fullMatrix<double> dzdX(nSamplingPnts, numEl);
643 fullMatrix<double> dxdY(nSamplingPnts, numEl);
644 fullMatrix<double> dydY(nSamplingPnts, numEl);
645 fullMatrix<double> dzdY(nSamplingPnts, numEl);
646 fullMatrix<double> dxdZ(nSamplingPnts, numEl);
647 fullMatrix<double> dydZ(nSamplingPnts, numEl);
648 fullMatrix<double> dzdZ(nSamplingPnts, numEl);
649 dSMat_dX.mult(nodesX, dxdX);
650 dSMat_dX.mult(nodesY, dydX);
651 dSMat_dX.mult(nodesZ, dzdX);
652 dSMat_dY.mult(nodesX, dxdY);
653 dSMat_dY.mult(nodesY, dydY);
654 dSMat_dY.mult(nodesZ, dzdY);
655 dSMat_dZ.mult(nodesX, dxdZ);
656 dSMat_dZ.mult(nodesY, dydZ);
657 dSMat_dZ.mult(nodesZ, dzdZ);
658 for(int iEl = 0; iEl < numEl; iEl++) {
659 for(int i = 0; i < nSamplingPnts; i++)
660 jacobian(i, iEl) = calcDet3D(dxdX(i, iEl), dxdY(i, iEl), dxdZ(i, iEl),
661 dydX(i, iEl), dydY(i, iEl), dydZ(i, iEl),
662 dzdX(i, iEl), dzdY(i, iEl), dzdZ(i, iEl));
663 if(scaling) {
664 fullMatrix<double> nodesXYZ(numPrimMapNodes, 3);
665 for(int i = 0; i < numPrimMapNodes; i++) {
666 nodesXYZ(i, 0) = nodesX(i, iEl);
667 nodesXYZ(i, 1) = nodesY(i, iEl);
668 nodesXYZ(i, 2) = nodesZ(i, iEl);
669 }
670 const double scale = 1. / getPrimJac3D(nodesXYZ);
671 for(int i = 0; i < nSamplingPnts; i++) jacobian(i, iEl) *= scale;
672 }
673 }
674 } break;
675 }
676 }
677
678 // Calculate the (signed) Jacobian determinant (in short, J) and its partial
679 // derivatives w.r.t. nodes coordinates for the element defined by
680 // the given node positions, with given normal vectors to straight element
681 // for regularization of 1D and 2D elements).
682 // Sampling points depend on the input matrices of shape function partial
683 // derivatives 'dSMat_d*', and only the 'nSamplingPnts' first of them
684 // are computed.
685 // The result is written in the matrix 'JDJ' which should be of size at
686 // least "nSamplingPnts x (3 * numMapNodes + 1)".
687 // For each sampling point, a row of 'JDJ' is filled with:
688 // - the partial derivatives of J w.r.t. the x component of the nodes
689 // - the partial derivatives of J w.r.t. the y component of the nodes
690 // - the partial derivatives of J w.r.t. the z component of the nodes
691 // - J
692 // NB: (x, y, z) are the physical coordinates and (X, Y, Z) are the reference
693 // coordinates
getSignedJacAndGradientsGeneral(int nSamplingPnts,const fullMatrix<double> & dSMat_dX,const fullMatrix<double> & dSMat_dY,const fullMatrix<double> & dSMat_dZ,const fullMatrix<double> & nodesXYZ,const fullMatrix<double> & normals,fullMatrix<double> & JDJ) const694 void JacobianBasis::getSignedJacAndGradientsGeneral(
695 int nSamplingPnts, const fullMatrix<double> &dSMat_dX,
696 const fullMatrix<double> &dSMat_dY, const fullMatrix<double> &dSMat_dZ,
697 const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
698 fullMatrix<double> &JDJ) const
699 {
700 switch(_dim) {
701 case 0: {
702 for(int i = 0; i < nSamplingPnts; i++) {
703 for(int j = 0; j < numMapNodes; j++) {
704 JDJ(i, j) = 0.;
705 JDJ(i, j + 1 * numMapNodes) = 0.;
706 JDJ(i, j + 2 * numMapNodes) = 0.;
707 }
708 JDJ(i, 3 * numMapNodes) = 1.;
709 }
710 } break;
711 case 1: {
712 fullMatrix<double> dxyzdX(nSamplingPnts, 3);
713 dSMat_dX.mult(nodesXYZ, dxyzdX);
714 for(int i = 0; i < nSamplingPnts; i++) {
715 calcJDJ1D(dxyzdX(i, 0), normals(0, 0), normals(1, 0), dxyzdX(i, 1),
716 normals(0, 1), normals(1, 1), dxyzdX(i, 2), normals(0, 2),
717 normals(1, 2), i, numMapNodes, dSMat_dX, JDJ);
718 }
719 } break;
720 case 2: {
721 fullMatrix<double> dxyzdX(nSamplingPnts, 3);
722 fullMatrix<double> dxyzdY(nSamplingPnts, 3);
723 dSMat_dX.mult(nodesXYZ, dxyzdX);
724 dSMat_dY.mult(nodesXYZ, dxyzdY);
725 for(int i = 0; i < nSamplingPnts; i++) {
726 calcJDJ2D(dxyzdX(i, 0), dxyzdY(i, 0), normals(0, 0), dxyzdX(i, 1),
727 dxyzdY(i, 1), normals(0, 1), dxyzdX(i, 2), dxyzdY(i, 2),
728 normals(0, 2), i, numMapNodes, dSMat_dX, dSMat_dY, JDJ);
729 }
730 } break;
731 case 3: {
732 fullMatrix<double> dxyzdX(nSamplingPnts, 3);
733 fullMatrix<double> dxyzdY(nSamplingPnts, 3);
734 fullMatrix<double> dxyzdZ(nSamplingPnts, 3);
735 dSMat_dX.mult(nodesXYZ, dxyzdX);
736 dSMat_dY.mult(nodesXYZ, dxyzdY);
737 dSMat_dZ.mult(nodesXYZ, dxyzdZ);
738 for(int i = 0; i < nSamplingPnts; i++) {
739 calcJDJ3D(dxyzdX(i, 0), dxyzdY(i, 0), dxyzdZ(i, 0), dxyzdX(i, 1),
740 dxyzdY(i, 1), dxyzdZ(i, 1), dxyzdX(i, 2), dxyzdY(i, 2),
741 dxyzdZ(i, 2), i, numMapNodes, dSMat_dX, dSMat_dY, dSMat_dZ,
742 JDJ);
743 }
744 } break;
745 }
746 }
747
getSignedIdealJacAndGradientsGeneral(int nSamplingPnts,const fullMatrix<double> & dSMat_dX,const fullMatrix<double> & dSMat_dY,const fullMatrix<double> & dSMat_dZ,const fullMatrix<double> & nodesXYZ,const fullMatrix<double> & normals,fullMatrix<double> & JDJ) const748 void JacobianBasis::getSignedIdealJacAndGradientsGeneral(
749 int nSamplingPnts, const fullMatrix<double> &dSMat_dX,
750 const fullMatrix<double> &dSMat_dY, const fullMatrix<double> &dSMat_dZ,
751 const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
752 fullMatrix<double> &JDJ) const
753 {
754 getSignedJacAndGradientsGeneral(nSamplingPnts, dSMat_dX, dSMat_dY, dSMat_dZ,
755 nodesXYZ, normals, JDJ);
756 }
757
getMetricMinAndGradients(const fullMatrix<double> & nodesXYZ,const fullMatrix<double> & nodesXYZStraight,fullVector<double> & lambdaJ,fullMatrix<double> & gradLambdaJ) const758 void JacobianBasis::getMetricMinAndGradients(
759 const fullMatrix<double> &nodesXYZ,
760 const fullMatrix<double> &nodesXYZStraight, fullVector<double> &lambdaJ,
761 fullMatrix<double> &gradLambdaJ) const
762 {
763 // jacobian of the straight elements (only triangles for now)
764 SPoint3 v0(nodesXYZ(0, 0), nodesXYZ(0, 1), nodesXYZ(0, 2));
765 SPoint3 v1(nodesXYZ(1, 0), nodesXYZ(1, 1), nodesXYZ(1, 2));
766 SPoint3 v2(nodesXYZ(2, 0), nodesXYZ(2, 1), nodesXYZ(2, 2));
767 SPoint3 *IXYZ[3] = {&v0, &v1, &v2};
768 double jaci[2][2] = {
769 {IXYZ[1]->x() - IXYZ[0]->x(), IXYZ[2]->x() - IXYZ[0]->x()},
770 {IXYZ[1]->y() - IXYZ[0]->y(), IXYZ[2]->y() - IXYZ[0]->y()}};
771 double invJaci[2][2];
772 inv2x2(jaci, invJaci);
773
774 for(int l = 0; l < numSamplingPnts; l++) {
775 double jac[2][2] = {{0., 0.}, {0., 0.}};
776 for(int i = 0; i < numMapNodes; i++) {
777 const double &dPhidX = _gradBasis->dShapeMat_dX(l, i);
778 const double &dPhidY = _gradBasis->dShapeMat_dY(l, i);
779 const double dpsidx = dPhidX * invJaci[0][0] + dPhidY * invJaci[1][0];
780 const double dpsidy = dPhidX * invJaci[0][1] + dPhidY * invJaci[1][1];
781 jac[0][0] += nodesXYZ(i, 0) * dpsidx;
782 jac[0][1] += nodesXYZ(i, 0) * dpsidy;
783 jac[1][0] += nodesXYZ(i, 1) * dpsidx;
784 jac[1][1] += nodesXYZ(i, 1) * dpsidy;
785 }
786 const double dxdx = jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1];
787 const double dydy = jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1];
788 const double dxdy = jac[0][0] * jac[1][0] + jac[0][1] * jac[1][1];
789 const double sqr = sqrt((dxdx - dydy) * (dxdx - dydy) + 4 * dxdy * dxdy);
790 const double osqr = sqr > 1e-8 ? 1 / sqr : 0;
791 lambdaJ(l) = 0.5 * (dxdx + dydy - sqr);
792 const double axx =
793 (1 - (dxdx - dydy) * osqr) * jac[0][0] - 2 * dxdy * osqr * jac[1][0];
794 const double axy =
795 (1 - (dxdx - dydy) * osqr) * jac[0][1] - 2 * dxdy * osqr * jac[1][1];
796 const double ayx =
797 -2 * dxdy * osqr * jac[0][0] + (1 - (dydy - dxdx) * osqr) * jac[1][0];
798 const double ayy =
799 -2 * dxdy * osqr * jac[0][1] + (1 - (dydy - dxdx) * osqr) * jac[1][1];
800 const double axixi = axx * invJaci[0][0] + axy * invJaci[0][1];
801 const double aetaeta = ayx * invJaci[1][0] + ayy * invJaci[1][1];
802 const double aetaxi = ayx * invJaci[0][0] + ayy * invJaci[0][1];
803 const double axieta = axx * invJaci[1][0] + axy * invJaci[1][1];
804 for(int i = 0; i < numMapNodes; i++) {
805 const double &dPhidX = _gradBasis->dShapeMat_dX(l, i);
806 const double &dPhidY = _gradBasis->dShapeMat_dY(l, i);
807 gradLambdaJ(l, i + 0 * numMapNodes) = axixi * dPhidX + axieta * dPhidY;
808 gradLambdaJ(l, i + 1 * numMapNodes) = aetaxi * dPhidX + aetaeta * dPhidY;
809 }
810 }
811 }
812
jacobianOrder(int tag)813 int JacobianBasis::jacobianOrder(int tag)
814 {
815 const int parentType = ElementType::getParentType(tag);
816 const int order = ElementType::getOrder(tag);
817 return jacobianOrder(parentType, order);
818 }
819
jacobianOrder(int parentType,int order)820 int JacobianBasis::jacobianOrder(int parentType, int order)
821 {
822 switch(parentType) {
823 case TYPE_PNT: return 0;
824 case TYPE_LIN: return order - 1;
825 case TYPE_TRI: return 2 * order - 2;
826 case TYPE_QUA: return 2 * order - 1;
827 case TYPE_TET: return 3 * order - 3;
828 case TYPE_PRI: return 3 * order - 1;
829 case TYPE_HEX: return 3 * order - 1;
830 case TYPE_PYR:
831 return 3 * order - 3;
832 // note : for the pyramid, the jacobian space is
833 // different from the space of the mapping
834 default:
835 Msg::Error("Unknown element type %d, return order 0", parentType);
836 return 0;
837 }
838 }
839
jacobianMatrixSpace(int type,int order)840 FuncSpaceData JacobianBasis::jacobianMatrixSpace(int type, int order)
841 {
842 if(type == TYPE_PYR) {
843 Msg::Error("jacobianMatrixSpace not yet implemented for pyramids");
844 return FuncSpaceData(type, false, 1, 0, false);
845 }
846 int jacOrder = -1;
847 switch(type) {
848 case TYPE_PNT: jacOrder = 0; break;
849 case TYPE_LIN:
850 case TYPE_TRI:
851 case TYPE_TET: jacOrder = order - 1; break;
852 case TYPE_QUA:
853 case TYPE_PRI:
854 case TYPE_HEX: jacOrder = order; break;
855 default:
856 Msg::Error("Unknown element type %d, return default space", type);
857 return FuncSpaceData();
858 }
859 return FuncSpaceData(type, jacOrder, false);
860 }
861