v0.9.0
UserDataOperators.cpp
Go to the documentation of this file.
1 /** \file ForcesAndSourcesCore.cpp
2 
3 \brief Implementation of Elements on Entities for Forces and Sources
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 namespace MoFEM {
22 
24 OpCalculateJacForFace::doWork(int side, EntityType type,
26 
28 
29  auto cal_jac_on_tri = [&]() {
31  if (type == MBVERTEX) {
32  VectorDouble &coords = getCoords();
33  double *coords_ptr = &*coords.data().begin();
34  double j00 = 0, j01 = 0, j10 = 0, j11 = 0;
35  // this is triangle, derivative of nodal shape functions is constant.
36  // So only need to do one node.
37  for (auto n : {0, 1, 2}) {
38  j00 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 0];
39  j01 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 1];
40  j10 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 0];
41  j11 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 1];
42  }
43  size_t nb_gauss_pts = getGaussPts().size2();
44  jac.resize(4, nb_gauss_pts, false);
46  &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0));
47  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
48  t_jac(0, 0) = j00;
49  t_jac(0, 1) = j01;
50  t_jac(1, 0) = j10;
51  t_jac(1, 1) = j11;
52  }
53  }
55  };
56 
57  auto cal_jac_on_quad = [&]() {
59  if (type == MBVERTEX) {
60  VectorDouble &coords = getCoords();
61  double *coords_ptr = &*coords.data().begin();
62  const size_t nb_integration_pts = getGaussPts().size2();
63  jac.resize(4, nb_integration_pts, false);
65  &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0));
66  double *ksi_ptr = &getGaussPts()(0, 0);
67  double *zeta_ptr = &getGaussPts()(1, 0);
68  for (size_t gg = 0; gg != nb_integration_pts;
69  ++gg, ++t_jac, ++ksi_ptr, ++zeta_ptr) {
70  const double &ksi = *ksi_ptr;
71  const double &zeta = *zeta_ptr;
72  jac(0, 0) = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0x(zeta) +
73  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1x(zeta) +
74  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2x(zeta) +
75  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3x(zeta);
76  jac(0, 1) = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0y(ksi) +
77  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1y(ksi) +
78  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2y(ksi) +
79  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3y(ksi);
80  jac(1, 0) = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0x(zeta) +
81  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1x(zeta) +
82  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2x(zeta) +
83  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3x(zeta);
84  jac(1, 1) = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0y(ksi) +
85  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1y(ksi) +
86  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2y(ksi) +
87  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3y(ksi);
88  }
89  }
91  };
92 
93  switch (getNumeredEntFiniteElementPtr()->getEntType()) {
94  case MBTRI:
95  CHKERR cal_jac_on_tri();
96  break;
97  case MBQUAD:
98  CHKERR cal_jac_on_quad();
99  break;
100  default:
101  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
102  "Operator not implemented for this entity type");
103  };
104 
105  doVertices = true;
106  doEdges = false;
107  doQuads = false;
108  doTris = false;
109  doTets = false;
110  doPrisms = false;
111 
113 }
114 
116 OpCalculateInvJacForFace::doWork(int side, EntityType type,
118 
120 
121  auto cal_inv_jac_on_tri = [&]() {
123  if (type == MBVERTEX) {
124  VectorDouble &coords = getCoords();
125  double *coords_ptr = &*coords.data().begin();
126  double j00 = 0, j01 = 0, j10 = 0, j11 = 0;
127 
128  // this is triangle, derivative of nodal shape functions is constant.
129  // So only need to do one node.
130  for (auto n : {0, 1, 2}) {
131  j00 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 0];
132  j01 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 1];
133  j10 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 0];
134  j11 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 1];
135  }
136  const double det = j00 * j11 - j01 * j10;
137 
138  size_t nb_gauss_pts = getGaussPts().size2();
139  invJac.resize(4, nb_gauss_pts, false);
141  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
142  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
143  t_inv_jac(0, 0) = j11 / det;
144  t_inv_jac(0, 1) = -j01 / det;
145  t_inv_jac(1, 0) = -j10 / det;
146  t_inv_jac(1, 1) = j00 / det;
147  }
148  }
150  };
151 
152  auto cal_inv_jac_on_quad = [&]() {
154  if (type == MBVERTEX) {
155  VectorDouble &coords = getCoords();
156  double *coords_ptr = &*coords.data().begin();
157  size_t nb_integration_pts = getGaussPts().size2();
158  invJac.resize(4, nb_integration_pts, false);
160  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
161  double *ksi_ptr = &getGaussPts()(0, 0);
162  double *zeta_ptr = &getGaussPts()(1, 0);
163  for (size_t gg = 0; gg != nb_integration_pts;
164  ++gg, ++t_inv_jac, ++ksi_ptr, ++zeta_ptr) {
165  const double &ksi = *ksi_ptr;
166  const double &zeta = *zeta_ptr;
167  double j00 = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0x(zeta) +
168  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1x(zeta) +
169  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2x(zeta) +
170  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3x(zeta);
171  double j01 = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0y(ksi) +
172  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1y(ksi) +
173  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2y(ksi) +
174  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3y(ksi);
175  double j10 = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0x(zeta) +
176  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1x(zeta) +
177  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2x(zeta) +
178  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3x(zeta);
179  double j11 = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0y(ksi) +
180  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1y(ksi) +
181  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2y(ksi) +
182  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3y(ksi);
183  double det = j00 * j11 - j01 * j10;
184  t_inv_jac(0, 0) = j11 / det;
185  t_inv_jac(0, 1) = -j01 / det;
186  t_inv_jac(1, 0) = -j10 / det;
187  t_inv_jac(1, 1) = j00 / det;
188  }
189  }
191  };
192 
193  switch (getNumeredEntFiniteElementPtr()->getEntType()) {
194  case MBTRI:
195  CHKERR cal_inv_jac_on_tri();
196  break;
197  case MBQUAD:
198  CHKERR cal_inv_jac_on_quad();
199  break;
200  default:
201  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
202  "Operator not implemented for this entity type");
203  };
204 
205  doVertices = true;
206  doEdges = false;
207  doQuads = false;
208  doTris = false;
209  doTets = false;
210  doPrisms = false;
211 
213 }
214 
216 OpSetInvJacH1ForFace::doWork(int side, EntityType type,
219 
220  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI &&
221  getNumeredEntFiniteElementPtr()->getEntType() != MBQUAD)
222  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
223  "This operator can be used only with element which is triangle");
224 
225  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
226 
227  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
228 
229  unsigned int nb_functions = data.getN(base).size2();
230  if (nb_functions) {
231  unsigned int nb_gauss_pts = data.getN(base).size1();
232  diffNinvJac.resize(nb_gauss_pts, 2 * nb_functions, false);
233 
234  if (type != MBVERTEX) {
235  if (nb_functions != data.getDiffN(base).size2() / 2) {
236  SETERRQ2(
237  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
238  "data inconsistency nb_functions != data.diffN.size2()/2 ( %u != "
239  "%u/2 )",
240  nb_functions, data.getDiffN(base).size2());
241  }
242  }
243 
244  switch (type) {
245  case MBVERTEX:
246  case MBEDGE:
247  case MBTRI:
248  case MBQUAD: {
253  &diffNinvJac(0, 0), &diffNinvJac(0, 1));
255  &data.getDiffN(base)(0, 0), &data.getDiffN(base)(0, 1));
257  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
258  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
259  for (unsigned int dd = 0; dd != nb_functions; ++dd) {
260  t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
261  ++t_diff_n;
262  ++t_diff_n_ref;
263  }
264  }
265  data.getDiffN(base).data().swap(diffNinvJac.data());
266  } break;
267  default:
268  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
269  }
270  }
271  }
272 
274 }
275 
277 OpSetInvJacHcurlFace::doWork(int side, EntityType type,
280 
281  if (type != MBEDGE && type != MBTRI)
283 
284  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI)
285  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
286  "This operator can be used only with element which is triangle");
287 
291 
292  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
293 
294  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
295 
296  const unsigned int nb_base_functions = data.getDiffN(base).size2() / 6;
297  if (nb_base_functions) {
298  const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
299 
300  diffHcurlInvJac.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
301 
302  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
303  double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
305  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
306 
307  &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
308 
309  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1]);
310 
312  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
313 
314  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
315  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
316  t_inv_diff_n(i, j) = t_diff_n(i, k) * t_inv_jac(k, j);
317  ++t_diff_n;
318  ++t_inv_diff_n;
319  }
320  }
321 
322  data.getDiffN(base).data().swap(diffHcurlInvJac.data());
323  }
324  }
325 
327 }
328 
330 OpMakeHdivFromHcurl::doWork(int side, EntityType type,
333 
334  if (type != MBEDGE && type != MBTRI)
336 
337  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI)
338  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
339  "This operator can be used only with element which is triangle");
340 
341  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
342 
343  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
344 
345  const size_t nb_base_functions = data.getN(base).size2() / 3;
346  if (nb_base_functions) {
347 
348  const size_t nb_gauss_pts = data.getN(base).size1();
349 
350  auto t_n = data.getFTensor1N<3>(base);
351  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
352  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
353  for (size_t bb = 0; bb != nb_base_functions; ++bb) {
354 
355  const double a = t_n(0);
356  t_n(0) = -t_n(1);
357  t_n(1) = a;
358 
359  for (auto n : {0, 1}) {
360  const double b = t_diff_n(0, n);
361  t_diff_n(0, n) = -t_diff_n(1, n);
362  t_diff_n(1, n) = b;
363  }
364 
365  ++t_n;
366  ++t_diff_n;
367  }
368  }
369  }
370  }
371 
373 }
374 
376  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
377 
379 
380  if (type != MBEDGE && type != MBTRI)
382 
386 
387  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
388 
389  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
390 
391  const size_t nb_base_functions = data.getN(base).size2() / 3;
392  if (nb_base_functions) {
393 
394  const size_t nb_gauss_pts = data.getN(base).size1();
395  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
396  if (data.getN(base).size2() > 0) {
397  auto t_n = data.getFTensor1N<3>(base);
398  double *t_transformed_n_ptr = &*piolaN.data().begin();
400  t_transformed_n_ptr, // HVEC0
401  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
403  &jAc(0, 0), &jAc(1, 0), &jAc(2, 0), &jAc(3, 0));
404  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
405  double det;
406  CHKERR determinantTensor2by2(t_jac, det);
407  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
408  t_transformed_n(i) = t_jac(i, k) * t_n(k) / det;
409  ++t_n;
410  ++t_transformed_n;
411  }
412  }
413  data.getN(base).data().swap(piolaN.data());
414  }
415 
416  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
417  if (data.getDiffN(base).size2() > 0) {
418  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
419  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
421  t_transformed_diff_n(t_transformed_diff_n_ptr,
422  &t_transformed_diff_n_ptr[HVEC0_1],
423 
424  &t_transformed_diff_n_ptr[HVEC1_0],
425  &t_transformed_diff_n_ptr[HVEC1_1],
426 
427  &t_transformed_diff_n_ptr[HVEC2_0],
428  &t_transformed_diff_n_ptr[HVEC2_1]);
430  &jAc(0, 0), &jAc(1, 0), &jAc(2, 0), &jAc(3, 0));
431  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
432  double det;
433  CHKERR determinantTensor2by2(t_jac, det);
434  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
435  t_transformed_diff_n(i, k) = t_jac(i, j) * t_diff_n(j, k) / det;
436  ++t_diff_n;
437  ++t_transformed_diff_n;
438  }
439  }
440  data.getDiffN(base).data().swap(piolaDiffN.data());
441  }
442  }
443  }
444 
446 }
447 
449  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
451 
452  if (type != MBEDGE)
454 
455  const auto &edge_direction = getDirection();
456 
458  FTensor::Tensor1<double, 3> t_m(-edge_direction[1], edge_direction[0],
459  edge_direction[2]);
460  const double l0 = t_m(i) * t_m(i);
461 
462  std::vector<double> l1;
463  {
464  int nb_gauss_pts = getTangetAtGaussPts().size1();
465  if (nb_gauss_pts) {
466  l1.resize(nb_gauss_pts);
467  const auto &edge_direction = getTangetAtGaussPts();
469  &edge_direction(0, 0), &edge_direction(0, 1), &edge_direction(0, 2));
470  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
471  l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
472  ++t_m_at_pts;
473  }
474  }
475  }
476 
477  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
478 
479  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
480  size_t nb_gauss_pts = data.getN(base).size1();
481  size_t nb_dofs = data.getN(base).size2() / 3;
482  if (nb_gauss_pts > 0 && nb_dofs > 0) {
483 
484  auto t_h_div = data.getFTensor1N<3>(base);
485 
486  size_t cc = 0;
487  const auto &edge_direction_at_gauss_pts = getTangetAtGaussPts();
488  if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
489 
491  &edge_direction_at_gauss_pts(0, 1),
492  &edge_direction_at_gauss_pts(0, 0),
493  &edge_direction_at_gauss_pts(0, 2));
494 
495  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
496  const double l0 = l1[gg];
497  for (int ll = 0; ll != nb_dofs; ++ll) {
498  const double val = t_h_div(0);
499  const double a = val / l0;
500  t_h_div(i) = t_m_at_pts(i) * a;
501  t_h_div(0) *= -1;
502  ++t_h_div;
503  ++cc;
504  }
505  ++t_m_at_pts;
506  }
507 
508  } else {
509 
510  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
511  for (int ll = 0; ll != nb_dofs; ll++) {
512  const double val = t_h_div(0);
513  const double a = val / l0;
514  t_h_div(i) = t_m(i) * a;
515  ++t_h_div;
516  ++cc;
517  }
518  }
519  }
520 
521  if (cc != nb_gauss_pts * nb_dofs)
522  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Data inconsistency");
523  }
524  }
525 
527 }
528 
530 OpCalculateInvJacForFatPrism::doWork(int side, EntityType type,
533 
534  if (type == MBVERTEX) {
535 
536  VectorDouble &coords = getCoords();
537  double *coords_ptr = &*coords.data().begin();
538 
539  const int nb_gauss_pts = data.getN(NOBASE).size1();
540  auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
541  invJac.resize(9, nb_gauss_pts, false);
542  invJac.clear();
543  auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
544 
548 
549  double &vol = getVolume();
550  vol = 0;
551 
552  for (int gg = 0; gg != nb_gauss_pts; gg++) {
553 
554  FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
555  &coords_ptr[2], 3);
556  t_jac(i, j) = 0;
557  for (int bb = 0; bb != 6; bb++) {
558  t_jac(i, j) += t_coords(i) * t_diff_n(j);
559  ++t_diff_n;
560  ++t_coords;
561  }
562 
563  double det;
564  CHKERR determinantTensor3by3(t_jac, det);
565  CHKERR invertTensor3by3(t_jac, det, t_inv_jac);
566  ++t_inv_jac;
567 
568  vol += 0.5 * det * getGaussPts()(3, gg);
569  }
570  }
571 
572  doVertices = true;
573  doEdges = false;
574  doQuads = false;
575  doTris = false;
576  doTets = false;
577  doPrisms = false;
578 
580 }
581 
583 OpSetInvJacH1ForFatPrism::doWork(int side, EntityType type,
586 
587  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
588 
589  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
590  if (data.getN(base).size2() == 0)
591  continue;
592 
593  const int nb_gauss_pts = data.getN(base).size1();
594  auto t_diff_n = data.getFTensor1DiffN<3>(base);
595  diffNinvJac.resize(data.getDiffN(base).size1(), data.getDiffN(base).size2(),
596  false);
597 
600 
602  &diffNinvJac(0, 0), &diffNinvJac(0, 1), &diffNinvJac(0, 2));
603  auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
604 
605  const int nb_dofs = data.getN(base).size2();
606  for (int gg = 0; gg != nb_gauss_pts; gg++) {
607  for (int bb = 0; bb != nb_dofs; bb++) {
608  t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
609  ++t_inv_diff_n;
610  ++t_diff_n;
611  }
612  ++t_inv_jac;
613  }
614 
615  data.getDiffN(base).data().swap(diffNinvJac.data());
616  }
617 
619 }
620 
622 OpCalculateInvJacForFlatPrism::doWork(int side, EntityType type,
624 
626 
627  if (type == MBVERTEX) {
628 
629  VectorDouble &coords = getCoords();
630  double *coords_ptr = &*coords.data().begin();
631  double diff_n[6];
632  CHKERR ShapeDiffMBTRI(diff_n);
633  double j00_f3, j01_f3, j10_f3, j11_f3;
634  for (int gg = 0; gg < 1; gg++) {
635  // this is triangle, derivative of nodal shape functions is constant.
636  // So only need to do one node.
637  j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
638  j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
639  j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
640  j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
641  }
642  double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
643  invJacF3.resize(2, 2, false);
644  invJacF3(0, 0) = j11_f3 / det_f3;
645  invJacF3(0, 1) = -j01_f3 / det_f3;
646  invJacF3(1, 0) = -j10_f3 / det_f3;
647  invJacF3(1, 1) = j00_f3 / det_f3;
648  }
649 
650  doVertices = true;
651  doEdges = false;
652  doQuads = false;
653  doTris = false;
654  doTets = false;
655  doPrisms = false;
656 
658 }
659 
661 OpSetInvJacH1ForFlatPrism::doWork(int side, EntityType type,
664 
665  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
666 
667  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
668 
669  unsigned int nb_dofs = data.getN(base).size2();
670  if (nb_dofs == 0)
672  unsigned int nb_gauss_pts = data.getN(base).size1();
673  diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs, false);
674 
675  if (type != MBVERTEX) {
676  if (nb_dofs != data.getDiffN(base).size2() / 2) {
677  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
678  "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
679  "%u/2 )",
680  nb_dofs, data.getDiffN(base).size2());
681  }
682  }
683 
684  switch (type) {
685  case MBVERTEX:
686  case MBEDGE:
687  case MBTRI: {
688  for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
689  for (unsigned int dd = 0; dd < nb_dofs; dd++) {
691  &*invJacF3.data().begin(), 2,
692  &data.getDiffN(base)(gg, 2 * dd), 1, 0,
693  &diffNinvJac(gg, 2 * dd), 1);
694  }
695  }
696  data.getDiffN(base).data().swap(diffNinvJac.data());
697  } break;
698  default:
699  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
700  }
701  }
702 
704 }
705 
706 } // namespace MoFEM
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
Definition: cblas_dgemv.c:11
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:415
user implemented approximation base
Definition: definitions.h:154
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:75
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
bool doVertices
If false skip vertices.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:69
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:70
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:74
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:428
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:396
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:73
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
FieldApproximationBase
approximation base
Definition: definitions.h:144
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
bool doEdges
If false skip edges.
MatrixDouble & getTangetAtGaussPts()
get tangent vector to edge curve at integration points
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:87
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:71
Data on single entity (This is passed as argument to DataOperator::doWork)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:76
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:72