gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
JacobianBasis.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include "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>
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
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).
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).
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
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.
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;
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 
182 }
183 
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 
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 
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 
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 
232  fullMatrix<double> &dSMat_dY,
233  fullMatrix<double> &dSMat_dZ)
234 {
235  calcMapFromIdealElement(type, dSMat_dX, dSMat_dY, dSMat_dZ);
236 }
237 
239  fullVector<double> &dSVec_dY,
240  fullVector<double> &dSVec_dZ)
241 {
242  calcMapFromIdealElement(type, dSVec_dX, dSVec_dY, dSVec_dZ);
243 }
244 
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 
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 
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 
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 
307 
308  delete[] barDPsi;
309 
310  // Compute "fast" Jacobian evaluation matrices (at 1st order nodes +
311  // barycenter)
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;
322  mapBasis->df(lagPointsFast, allDPsiFast);
323  numMapNodes = mapBasis->getNumShapeFunctions();
324 
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)
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)
386  fullMatrix<double> &result,
387  bool ideal) const
388 {
389  const fullVector<double> &gSX =
391  const fullVector<double> &gSY =
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
419  bool ideal) const
420 {
421  const fullVector<double> &gSX =
423  const fullVector<double> &gSY =
425  const fullVector<double> &gSZ =
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.
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
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
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 
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 
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 
814 {
815  const int parentType = ElementType::getParentType(tag);
816  const int order = ElementType::getOrder(tag);
817  return jacobianOrder(parentType, order);
818 }
819 
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 
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 }
GradientBasis::getIdealGradientsFromNodes
void getIdealGradientsFromNodes(const fullMatrix< double > &nodesCoord, fullMatrix< double > *dxyzdX, fullMatrix< double > *dxyzdY, fullMatrix< double > *dxyzdZ) const
Definition: JacobianBasis.cpp:184
JacobianBasis::matrixPrimJac2Jac
fullMatrix< double > matrixPrimJac2Jac
Definition: JacobianBasis.h:71
GradientBasis::_data
const FuncSpaceData _data
Definition: JacobianBasis.h:19
JacobianBasis::getMetricMinAndGradients
void getMetricMinAndGradients(const fullMatrix< double > &nodesXYZ, const fullMatrix< double > &nodesXYZStraight, fullVector< double > &lambdaJ, fullMatrix< double > &gradLambdaJ) const
Definition: JacobianBasis.cpp:758
sqr
static int sqr(int x)
Definition: gl2gif.cpp:266
JacobianBasis::dFastShapeMat_dY
fullMatrix< double > dFastShapeMat_dY
Definition: JacobianBasis.h:66
JacobianBasis::numMapNodes
int numMapNodes
Definition: JacobianBasis.h:73
nodalBasis::points
fullMatrix< double > points
Definition: nodalBasis.h:16
JacobianBasis::JacobianBasis
JacobianBasis(int elementTag, FuncSpaceData)
Definition: JacobianBasis.cpp:252
fullVector< double >
TYPE_LIN
#define TYPE_LIN
Definition: GmshDefines.h:65
JacobianBasis::_dim
const int _dim
Definition: JacobianBasis.h:65
nodalBasis::df
virtual void df(double u, double v, double w, double grads[][3]) const =0
gmshGenerateOrderedPoints
void gmshGenerateOrderedPoints(FuncSpaceData data, fullMatrix< double > &points, bool onBezierDomain)
Definition: pointsGenerators.cpp:1080
nodalBasis.h
fullMatrix::setAll
void setAll(const scalar &m)
Definition: fullMatrix.h:422
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
GradientBasis::dShapeMat_dY
fullMatrix< double > dShapeMat_dY
Definition: JacobianBasis.h:14
SPoint3
Definition: SPoint3.h:14
GradientBasis::getAllGradientsFromNodes
void getAllGradientsFromNodes(const fullMatrix< double > &nodesCoord, fullMatrix< double > &dxyzdXYZ) const
Definition: JacobianBasis.cpp:203
TYPE_PNT
#define TYPE_PNT
Definition: GmshDefines.h:64
JacobianBasis::dPrimBaryIdealShape_dY
fullVector< double > dPrimBaryIdealShape_dY
Definition: JacobianBasis.h:68
JacobianBasis::_gradBasis
const GradientBasis * _gradBasis
Definition: JacobianBasis.h:62
fullVector::scale
void scale(const scalar s)
Definition: fullMatrix.h:144
TYPE_TRI
#define TYPE_TRI
Definition: GmshDefines.h:66
fullVector::resize
bool resize(int r, bool resetValue=true)
Definition: fullMatrix.h:103
GradientBasis::mapFromIdealElement
void mapFromIdealElement(fullMatrix< double > &dxyzdX, fullMatrix< double > &dxyzdY, fullMatrix< double > &dxyzdZ) const
Definition: JacobianBasis.h:39
JacobianBasis::getPrimNormals1D
double getPrimNormals1D(const fullMatrix< double > &nodesXYZ, fullMatrix< double > &result) const
Definition: JacobianBasis.cpp:339
JacobianBasis::_elementTag
const int _elementTag
Definition: JacobianBasis.h:63
TYPE_PRI
#define TYPE_PRI
Definition: GmshDefines.h:70
GradientBasis::dShapeMat_dX
fullMatrix< double > dShapeMat_dX
Definition: JacobianBasis.h:14
cTet
static const double cTet
Definition: qualityMeasuresJacobian.cpp:20
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
fullMatrix< double >
GradientBasis::dShapeMat_dZ
fullMatrix< double > dShapeMat_dZ
Definition: JacobianBasis.h:14
BasisFactory::getNodalBasis
static const nodalBasis * getNodalBasis(int tag)
Definition: BasisFactory.cpp:23
nodalBasis::f
virtual void f(double u, double v, double w, double *sf) const =0
JacobianBasis::numSamplingPntsFast
int numSamplingPntsFast
Definition: JacobianBasis.h:74
JacobianBasis::jacobianOrder
static int jacobianOrder(int tag)
Definition: JacobianBasis.cpp:813
JacobianBasis::dPrimBaryShape_dY
fullVector< double > dPrimBaryShape_dY
Definition: JacobianBasis.h:67
JacobianBasis::getSignedIdealJacAndGradientsGeneral
void 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) const
Definition: JacobianBasis.cpp:748
JacobianBasis::dPrimBaryIdealShape_dX
fullVector< double > dPrimBaryIdealShape_dX
Definition: JacobianBasis.h:68
ElementType::getType
int getType(int parentType, int order, bool serendip=false)
Definition: ElementType.cpp:757
JacobianBasis::dFastShapeMat_dX
fullMatrix< double > dFastShapeMat_dX
Definition: JacobianBasis.h:66
GradientBasis::getAllIdealGradientsFromNodes
void getAllIdealGradientsFromNodes(const fullMatrix< double > &nodesCoord, fullMatrix< double > &dxyzdXYZ) const
Definition: JacobianBasis.cpp:217
JacobianBasis::numPrimMapNodes
int numPrimMapNodes
Definition: JacobianBasis.h:73
Numeric.h
GradientBasis::getGradientsFromNodes
void getGradientsFromNodes(const fullMatrix< double > &nodesCoord, fullMatrix< double > *dxyzdX, fullMatrix< double > *dxyzdY, fullMatrix< double > *dxyzdZ) const
Definition: JacobianBasis.cpp:193
BasisFactory::getGradientBasis
static const GradientBasis * getGradientBasis(int tag, FuncSpaceData)
Definition: BasisFactory.cpp:105
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
gmshGeneratePoints
void gmshGeneratePoints(FuncSpaceData data, fullMatrix< double > &points)
Definition: pointsGenerators.cpp:18
JacobianBasis.h
JacobianBasis::numSamplingPnts
int numSamplingPnts
Definition: JacobianBasis.h:72
GradientBasis::GradientBasis
GradientBasis(int elementTag, FuncSpaceData)
Definition: JacobianBasis.cpp:146
FuncSpaceData::getType
int getType() const
Definition: FuncSpaceData.h:65
cPyr
static const double cPyr
Definition: qualityMeasuresJacobian.cpp:21
GradientBasis::dShapeIdealMat_dZ
fullMatrix< double > dShapeIdealMat_dZ
Definition: JacobianBasis.h:15
cTri
static const double cTri
Definition: qualityMeasuresJacobian.cpp:19
fullMatrix::size2
int size2() const
Definition: fullMatrix.h:275
JacobianBasis::dFastShapeMat_dZ
fullMatrix< double > dFastShapeMat_dZ
Definition: JacobianBasis.h:66
JacobianBasis::getPrimNormal2D
double getPrimNormal2D(const fullMatrix< double > &nodesXYZ, fullMatrix< double > &result, bool ideal=false) const
Definition: JacobianBasis.cpp:385
TYPE_PYR
#define TYPE_PYR
Definition: GmshDefines.h:69
FuncSpaceData
Definition: FuncSpaceData.h:16
JacobianBasis::dPrimBaryShape_dZ
fullVector< double > dPrimBaryShape_dZ
Definition: JacobianBasis.h:67
inv2x2
double inv2x2(double mat[2][2], double inv[2][2])
Definition: Numeric.cpp:214
TYPE_QUA
#define TYPE_QUA
Definition: GmshDefines.h:67
GradientBasis::_elementTag
const int _elementTag
Definition: JacobianBasis.h:18
fullMatrix::setAsProxy
void setAsProxy(const fullMatrix< scalar > &original)
Definition: fullMatrix.h:335
nodalBasis
Definition: nodalBasis.h:12
JacobianBasis::jacobianMatrixSpace
static FuncSpaceData jacobianMatrixSpace(int type, int order)
Definition: JacobianBasis.cpp:840
ElementType::getDimension
int getDimension(int type)
Definition: ElementType.cpp:297
fullMatrix::size1
int size1() const
Definition: fullMatrix.h:274
ElementType::getOrder
int getOrder(int type)
Definition: ElementType.cpp:158
GradientBasis::dShapeIdealMat_dY
fullMatrix< double > dShapeIdealMat_dY
Definition: JacobianBasis.h:15
JacobianBasis::getSignedJacAndGradientsGeneral
void 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) const
Definition: JacobianBasis.cpp:694
TYPE_HEX
#define TYPE_HEX
Definition: GmshDefines.h:71
GradientBasis::dShapeIdealMat_dX
fullMatrix< double > dShapeIdealMat_dX
Definition: JacobianBasis.h:15
fullMatrix::copy
void copy(const fullMatrix< scalar > &a, int i0, int ni, int j0, int nj, int desti0, int destj0)
Definition: fullMatrix.h:394
TYPE_TET
#define TYPE_TET
Definition: GmshDefines.h:68
JacobianBasis::dPrimBaryShape_dX
fullVector< double > dPrimBaryShape_dX
Definition: JacobianBasis.h:67
JacobianBasis::getPrimJac3D
double getPrimJac3D(const fullMatrix< double > &nodesXYZ, bool ideal=false) const
Definition: JacobianBasis.cpp:418
fullMatrix::mult
void mult(const fullVector< scalar > &x, fullVector< scalar > &y) const
Definition: fullMatrix.h:487
JacobianBasis::dPrimBaryIdealShape_dZ
fullVector< double > dPrimBaryIdealShape_dZ
Definition: JacobianBasis.h:69
JacobianBasis::getJacobianGeneral
void 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 > *normals) const
Definition: JacobianBasis.cpp:448
pointsGenerators.h
JacobianBasis::numPrimSamplingPnts
int numPrimSamplingPnts
Definition: JacobianBasis.h:72
fullMatrix::resize
bool resize(int r, int c, bool resetValue=true)
Definition: fullMatrix.h:307
ElementType::getParentType
int getParentType(int type)
Definition: ElementType.cpp:10
nodalBasis::getNumShapeFunctions
virtual int getNumShapeFunctions() const =0
BasisFactory.h
scale
static void scale(std::vector< double > &x, double s)
Definition: ConjugateGradients.cpp:21