v0.10.0
UserDataOperators.cpp
Go to the documentation of this file.
1 /** \file UserDataOperators.cpp
2 
3 \brief Generic user data operators for evaluate fields, and other common
4 purposes.
5 
6 */
7 
8 /* This file is part of MoFEM.
9  * MoFEM is free software: you can redistribute it and/or modify it under
10  * the terms of the GNU Lesser General Public License as published by the
11  * Free Software Foundation, either version 3 of the License, or (at your
12  * option) any later version.
13  *
14  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17  * License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21 
22 namespace MoFEM {
23 
25 OpCalculateJacForFace::doWork(int side, EntityType type,
27 
29 
30  auto cal_jac_on_tri = [&]() {
32  if (type == MBVERTEX) {
33  VectorDouble &coords = getCoords();
34  double *coords_ptr = &*coords.data().begin();
35  double j00 = 0, j01 = 0, j10 = 0, j11 = 0;
36  // this is triangle, derivative of nodal shape functions is constant.
37  // So only need to do one node.
38  for (auto n : {0, 1, 2}) {
39  j00 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 0];
40  j01 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 1];
41  j10 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 0];
42  j11 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 1];
43  }
44  size_t nb_gauss_pts = getGaussPts().size2();
45  jac.resize(4, nb_gauss_pts, false);
47  &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0));
48  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
49  t_jac(0, 0) = j00;
50  t_jac(0, 1) = j01;
51  t_jac(1, 0) = j10;
52  t_jac(1, 1) = j11;
53  }
54  }
56  };
57 
58  auto cal_jac_on_quad = [&]() {
60  if (type == MBVERTEX) {
61  VectorDouble &coords = getCoords();
62  double *coords_ptr = &*coords.data().begin();
63  const size_t nb_integration_pts = getGaussPts().size2();
64  jac.resize(4, nb_integration_pts, false);
66  &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0));
67  double *ksi_ptr = &getGaussPts()(0, 0);
68  double *zeta_ptr = &getGaussPts()(1, 0);
69  for (size_t gg = 0; gg != nb_integration_pts;
70  ++gg, ++t_jac, ++ksi_ptr, ++zeta_ptr) {
71  const double &ksi = *ksi_ptr;
72  const double &zeta = *zeta_ptr;
73  jac(0, 0) = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0x(zeta) +
74  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1x(zeta) +
75  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2x(zeta) +
76  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3x(zeta);
77  jac(0, 1) = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0y(ksi) +
78  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1y(ksi) +
79  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2y(ksi) +
80  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3y(ksi);
81  jac(1, 0) = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0x(zeta) +
82  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1x(zeta) +
83  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2x(zeta) +
84  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3x(zeta);
85  jac(1, 1) = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0y(ksi) +
86  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1y(ksi) +
87  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2y(ksi) +
88  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3y(ksi);
89  }
90  }
92  };
93 
94  switch (getNumeredEntFiniteElementPtr()->getEntType()) {
95  case MBTRI:
96  CHKERR cal_jac_on_tri();
97  break;
98  case MBQUAD:
99  CHKERR cal_jac_on_quad();
100  break;
101  default:
102  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
103  "Operator not implemented for this entity type");
104  };
105 
106  doEntities[MBVERTEX] = true;
107  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
108 
110 }
111 
115 
117 
118  auto cal_inv_jac_on_tri = [&]() {
120  if (type == MBVERTEX) {
121  VectorDouble &coords = getCoords();
122  double *coords_ptr = &*coords.data().begin();
123  double j00 = 0, j01 = 0, j10 = 0, j11 = 0;
124 
125  // this is triangle, derivative of nodal shape functions is constant.
126  // So only need to do one node.
127  for (auto n : {0, 1, 2}) {
128  j00 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 0];
129  j01 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 1];
130  j10 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 0];
131  j11 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 1];
132  }
133  const double det = j00 * j11 - j01 * j10;
134 
135  size_t nb_gauss_pts = getGaussPts().size2();
136  invJac.resize(4, nb_gauss_pts, false);
138  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
139  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
140  t_inv_jac(0, 0) = j11 / det;
141  t_inv_jac(0, 1) = -j01 / det;
142  t_inv_jac(1, 0) = -j10 / det;
143  t_inv_jac(1, 1) = j00 / det;
144  }
145  }
147  };
148 
149  auto cal_inv_jac_on_quad = [&]() {
151  if (type == MBVERTEX) {
152  VectorDouble &coords = getCoords();
153  double *coords_ptr = &*coords.data().begin();
154  size_t nb_integration_pts = getGaussPts().size2();
155  invJac.resize(4, nb_integration_pts, false);
157  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
158  double *ksi_ptr = &getGaussPts()(0, 0);
159  double *zeta_ptr = &getGaussPts()(1, 0);
160  for (size_t gg = 0; gg != nb_integration_pts;
161  ++gg, ++t_inv_jac, ++ksi_ptr, ++zeta_ptr) {
162  const double &ksi = *ksi_ptr;
163  const double &zeta = *zeta_ptr;
164  double j00 = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0x(zeta) +
165  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1x(zeta) +
166  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2x(zeta) +
167  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3x(zeta);
168  double j01 = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0y(ksi) +
169  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1y(ksi) +
170  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2y(ksi) +
171  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3y(ksi);
172  double j10 = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0x(zeta) +
173  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1x(zeta) +
174  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2x(zeta) +
175  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3x(zeta);
176  double j11 = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0y(ksi) +
177  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1y(ksi) +
178  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2y(ksi) +
179  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3y(ksi);
180  double det = j00 * j11 - j01 * j10;
181  t_inv_jac(0, 0) = j11 / det;
182  t_inv_jac(0, 1) = -j01 / det;
183  t_inv_jac(1, 0) = -j10 / det;
184  t_inv_jac(1, 1) = j00 / det;
185  }
186  }
188  };
189 
190  switch (getNumeredEntFiniteElementPtr()->getEntType()) {
191  case MBTRI:
192  CHKERR cal_inv_jac_on_tri();
193  break;
194  case MBQUAD:
195  CHKERR cal_inv_jac_on_quad();
196  break;
197  default:
198  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
199  "Operator not implemented for this entity type");
200  };
201 
202  doEntities[MBVERTEX] = true;
203  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
204 
206 }
207 
209 OpSetInvJacH1ForFace::doWork(int side, EntityType type,
212 
213  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI &&
214  getNumeredEntFiniteElementPtr()->getEntType() != MBQUAD)
215  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
216  "This operator can be used only with element which is triangle");
217 
218  auto apply_transform = [&](MatrixDouble &diff_n) {
220  size_t nb_functions = diff_n.size2() / 2;
221  if (nb_functions) {
222  size_t nb_gauss_pts = diff_n.size1();
223  diffNinvJac.resize(nb_gauss_pts, 2 * nb_functions, false);
224 
225  switch (type) {
226  case MBVERTEX:
227  case MBEDGE:
228  case MBTRI:
229  case MBQUAD: {
234  &diffNinvJac(0, 0), &diffNinvJac(0, 1));
236  &diff_n(0, 0), &diff_n(0, 1));
238  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
239  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
240  for (size_t dd = 0; dd != nb_functions; ++dd) {
241  t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
242  ++t_diff_n;
243  ++t_diff_n_ref;
244  }
245  }
246  diff_n.data().swap(diffNinvJac.data());
247  } break;
248  default:
249  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
250  }
251  }
253  };
254 
255  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
256  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
257  CHKERR apply_transform(data.getDiffN(base));
258  }
259 
260  switch (type) {
261  case MBVERTEX:
262  for (auto &m : data.getBBDiffNMap())
263  CHKERR apply_transform(*(m.second));
264  break;
265  default:
266  for (auto &ptr : data.getBBDiffNByOrderArray())
267  if (ptr)
268  CHKERR apply_transform(*ptr);
269  }
270 
272 }
273 
275 OpSetInvJacHcurlFace::doWork(int side, EntityType type,
278 
279  if (type != MBEDGE && type != MBTRI)
281 
282  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI)
283  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
284  "This operator can be used only with element which is triangle");
285 
289 
290  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
291 
292  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
293 
294  const unsigned int nb_base_functions = data.getDiffN(base).size2() / 6;
295  if (nb_base_functions) {
296  const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
297 
298  diffHcurlInvJac.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
299 
300  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
301  double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
303  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
304 
305  &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
306 
307  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1]);
308 
310  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
311 
312  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
313  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
314  t_inv_diff_n(i, j) = t_diff_n(i, k) * t_inv_jac(k, j);
315  ++t_diff_n;
316  ++t_inv_diff_n;
317  }
318  }
319 
320  data.getDiffN(base).data().swap(diffHcurlInvJac.data());
321  }
322  }
323 
325 }
326 
328 OpMakeHdivFromHcurl::doWork(int side, EntityType type,
331 
332  if (type != MBEDGE && type != MBTRI)
334 
335  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI)
336  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
337  "This operator can be used only with element which is triangle");
338 
339  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
340 
341  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
342 
343  const size_t nb_base_functions = data.getN(base).size2() / 3;
344  if (nb_base_functions) {
345 
346  const size_t nb_gauss_pts = data.getN(base).size1();
347 
348  auto t_n = data.getFTensor1N<3>(base);
349  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
350  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
351  for (size_t bb = 0; bb != nb_base_functions; ++bb) {
352 
353  const double a = t_n(0);
354  t_n(0) = -t_n(1);
355  t_n(1) = a;
356 
357  for (auto n : {0, 1}) {
358  const double b = t_diff_n(0, n);
359  t_diff_n(0, n) = -t_diff_n(1, n);
360  t_diff_n(1, n) = b;
361  }
362 
363  ++t_n;
364  ++t_diff_n;
365  }
366  }
367  }
368  }
369 
371 }
372 
374  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
375 
377 
378  if (type != MBEDGE && type != MBTRI)
380 
384 
385  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
386 
387  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
388 
389  const size_t nb_base_functions = data.getN(base).size2() / 3;
390  if (nb_base_functions) {
391 
392  const size_t nb_gauss_pts = data.getN(base).size1();
393  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
394  if (data.getN(base).size2() > 0) {
395  auto t_n = data.getFTensor1N<3>(base);
396  double *t_transformed_n_ptr = &*piolaN.data().begin();
398  t_transformed_n_ptr, // HVEC0
399  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
401  &jAc(0, 0), &jAc(1, 0), &jAc(2, 0), &jAc(3, 0));
402  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
403  double det;
404  CHKERR determinantTensor2by2(t_jac, det);
405  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
406  t_transformed_n(i) = t_jac(i, k) * t_n(k) / det;
407  ++t_n;
408  ++t_transformed_n;
409  }
410  }
411  data.getN(base).data().swap(piolaN.data());
412  }
413 
414  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
415  if (data.getDiffN(base).size2() > 0) {
416  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
417  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
419  t_transformed_diff_n(t_transformed_diff_n_ptr,
420  &t_transformed_diff_n_ptr[HVEC0_1],
421 
422  &t_transformed_diff_n_ptr[HVEC1_0],
423  &t_transformed_diff_n_ptr[HVEC1_1],
424 
425  &t_transformed_diff_n_ptr[HVEC2_0],
426  &t_transformed_diff_n_ptr[HVEC2_1]);
428  &jAc(0, 0), &jAc(1, 0), &jAc(2, 0), &jAc(3, 0));
429  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
430  double det;
431  CHKERR determinantTensor2by2(t_jac, det);
432  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
433  t_transformed_diff_n(i, k) = t_jac(i, j) * t_diff_n(j, k) / det;
434  ++t_diff_n;
435  ++t_transformed_diff_n;
436  }
437  }
438  data.getDiffN(base).data().swap(piolaDiffN.data());
439  }
440  }
441  }
442 
444 }
445 
447  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
449 
450  if (type != MBEDGE)
452 
453  const auto &edge_direction = getDirection();
454 
456  FTensor::Tensor1<double, 3> t_m(-edge_direction[1], edge_direction[0],
457  edge_direction[2]);
458  const double l0 = t_m(i) * t_m(i);
459 
460  std::vector<double> l1;
461  {
462  int nb_gauss_pts = getTangetAtGaussPts().size1();
463  if (nb_gauss_pts) {
464  l1.resize(nb_gauss_pts);
465  const auto &edge_direction = getTangetAtGaussPts();
467  &edge_direction(0, 0), &edge_direction(0, 1), &edge_direction(0, 2));
468  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
469  l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
470  ++t_m_at_pts;
471  }
472  }
473  }
474 
475  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
476 
477  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
478  size_t nb_gauss_pts = data.getN(base).size1();
479  size_t nb_dofs = data.getN(base).size2() / 3;
480  if (nb_gauss_pts > 0 && nb_dofs > 0) {
481 
482  auto t_h_div = data.getFTensor1N<3>(base);
483 
484  size_t cc = 0;
485  const auto &edge_direction_at_gauss_pts = getTangetAtGaussPts();
486  if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
487 
489  &edge_direction_at_gauss_pts(0, 1),
490  &edge_direction_at_gauss_pts(0, 0),
491  &edge_direction_at_gauss_pts(0, 2));
492 
493  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
494  const double l0 = l1[gg];
495  for (int ll = 0; ll != nb_dofs; ++ll) {
496  const double val = t_h_div(0);
497  const double a = val / l0;
498  t_h_div(i) = t_m_at_pts(i) * a;
499  t_h_div(0) *= -1;
500  ++t_h_div;
501  ++cc;
502  }
503  ++t_m_at_pts;
504  }
505 
506  } else {
507 
508  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
509  for (int ll = 0; ll != nb_dofs; ll++) {
510  const double val = t_h_div(0);
511  const double a = val / l0;
512  t_h_div(i) = t_m(i) * a;
513  ++t_h_div;
514  ++cc;
515  }
516  }
517  }
518 
519  if (cc != nb_gauss_pts * nb_dofs)
520  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Data inconsistency");
521  }
522  }
523 
525 }
526 
528  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
530 
531  if (type == MBVERTEX) {
532 
533  VectorDouble &coords = getCoords();
534  double *coords_ptr = &*coords.data().begin();
535 
536  const int nb_gauss_pts = data.getN(NOBASE).size1();
537  auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
538 
542 
543  auto t_w = getFTensor0IntegrationWeight();
544  for (int gg = 0; gg != nb_gauss_pts; gg++) {
545 
546  FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
547  &coords_ptr[2], 3);
548  t_jac(i, j) = 0;
549  for (int bb = 0; bb != 6; bb++) {
550  t_jac(i, j) += t_coords(i) * t_diff_n(j);
551  ++t_diff_n;
552  ++t_coords;
553  }
554 
555  double det;
556  CHKERR determinantTensor3by3(t_jac, det);
557  t_w *= det / 2.;
558 
559  ++t_w;
560  }
561 
562  double &vol = getVolume();
563  auto t_w_scaled = getFTensor0IntegrationWeight();
564  for (int gg = 0; gg != nb_gauss_pts; gg++) {
565  t_w_scaled /= vol;
566  ++t_w_scaled;
567  }
568  }
569 
570  doEntities[MBVERTEX] = true;
571  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
572 
574 }
575 
580 
581  if (type == MBVERTEX) {
582 
583  VectorDouble &coords = getCoords();
584  double *coords_ptr = &*coords.data().begin();
585 
586  const int nb_gauss_pts = data.getN(NOBASE).size1();
587  auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
588  invJac.resize(9, nb_gauss_pts, false);
589  invJac.clear();
590  auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
591 
595 
596  auto t_w = getFTensor0IntegrationWeight();
597  for (int gg = 0; gg != nb_gauss_pts; gg++) {
598 
599  FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
600  &coords_ptr[2], 3);
601  t_jac(i, j) = 0;
602  for (int bb = 0; bb != 6; bb++) {
603  t_jac(i, j) += t_coords(i) * t_diff_n(j);
604  ++t_diff_n;
605  ++t_coords;
606  }
607 
608  double det;
609  CHKERR determinantTensor3by3(t_jac, det);
610  CHKERR invertTensor3by3(t_jac, det, t_inv_jac);
611  ++t_inv_jac;
612  }
613  }
614 
615  doEntities[MBVERTEX] = true;
616  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
617 
619 }
620 
625 
626  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
627 
628  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
629  if (data.getN(base).size2() == 0)
630  continue;
631 
632  const int nb_gauss_pts = data.getN(base).size1();
633  auto t_diff_n = data.getFTensor1DiffN<3>(base);
634  diffNinvJac.resize(data.getDiffN(base).size1(), data.getDiffN(base).size2(),
635  false);
636 
639 
641  &diffNinvJac(0, 0), &diffNinvJac(0, 1), &diffNinvJac(0, 2));
642  auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
643 
644  const int nb_dofs = data.getN(base).size2();
645  for (int gg = 0; gg != nb_gauss_pts; gg++) {
646  for (int bb = 0; bb != nb_dofs; bb++) {
647  t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
648  ++t_inv_diff_n;
649  ++t_diff_n;
650  }
651  ++t_inv_jac;
652  }
653 
654  data.getDiffN(base).data().swap(diffNinvJac.data());
655  }
656 
658 }
659 
663 
665 
666  if (type == MBVERTEX) {
667 
668  VectorDouble &coords = getCoords();
669  double *coords_ptr = &*coords.data().begin();
670  double diff_n[6];
671  CHKERR ShapeDiffMBTRI(diff_n);
672  double j00_f3, j01_f3, j10_f3, j11_f3;
673  for (int gg = 0; gg < 1; gg++) {
674  // this is triangle, derivative of nodal shape functions is constant.
675  // So only need to do one node.
676  j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
677  j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
678  j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
679  j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
680  }
681  double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
682  invJacF3.resize(2, 2, false);
683  invJacF3(0, 0) = j11_f3 / det_f3;
684  invJacF3(0, 1) = -j01_f3 / det_f3;
685  invJacF3(1, 0) = -j10_f3 / det_f3;
686  invJacF3(1, 1) = j00_f3 / det_f3;
687  }
688 
689  doEntities[MBVERTEX] = true;
690  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
691 
693 }
694 
699 
700  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
701 
702  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
703 
704  unsigned int nb_dofs = data.getN(base).size2();
705  if (nb_dofs == 0)
707  unsigned int nb_gauss_pts = data.getN(base).size1();
708  diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs, false);
709 
710  if (type != MBVERTEX) {
711  if (nb_dofs != data.getDiffN(base).size2() / 2) {
712  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
713  "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
714  "%u/2 )",
715  nb_dofs, data.getDiffN(base).size2());
716  }
717  }
718 
719  switch (type) {
720  case MBVERTEX:
721  case MBEDGE:
722  case MBTRI: {
723  for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
724  for (unsigned int dd = 0; dd < nb_dofs; dd++) {
726  &*invJacF3.data().begin(), 2,
727  &data.getDiffN(base)(gg, 2 * dd), 1, 0,
728  &diffNinvJac(gg, 2 * dd), 1);
729  }
730  }
731  data.getDiffN(base).data().swap(diffNinvJac.data());
732  } break;
733  default:
734  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
735  }
736  }
737 
739 }
740 
741 } // 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:571
user implemented approximation base
Definition: definitions.h:160
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
#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:509
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
#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.
T data[Tensor_Dim]
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
#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
static Index< 'n', 3 > n
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:584
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:552
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
static Index< 'm', 3 > m
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
static Index< 'i', 3 > i
FieldApproximationBase
approximation base
Definition: definitions.h:150
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
auto getFTensor0IntegrationWeight()
Get integration weights.
static Index< 'j', 3 > j
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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:604
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.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:71
static Index< 'k', 3 > k
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:415
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
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