v0.9.2
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  doEntities[MBVERTEX] = true;
106  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
107 
109 }
110 
114 
116 
117  auto cal_inv_jac_on_tri = [&]() {
119  if (type == MBVERTEX) {
120  VectorDouble &coords = getCoords();
121  double *coords_ptr = &*coords.data().begin();
122  double j00 = 0, j01 = 0, j10 = 0, j11 = 0;
123 
124  // this is triangle, derivative of nodal shape functions is constant.
125  // So only need to do one node.
126  for (auto n : {0, 1, 2}) {
127  j00 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 0];
128  j01 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 1];
129  j10 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 0];
130  j11 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 1];
131  }
132  const double det = j00 * j11 - j01 * j10;
133 
134  size_t nb_gauss_pts = getGaussPts().size2();
135  invJac.resize(4, nb_gauss_pts, false);
137  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
138  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
139  t_inv_jac(0, 0) = j11 / det;
140  t_inv_jac(0, 1) = -j01 / det;
141  t_inv_jac(1, 0) = -j10 / det;
142  t_inv_jac(1, 1) = j00 / det;
143  }
144  }
146  };
147 
148  auto cal_inv_jac_on_quad = [&]() {
150  if (type == MBVERTEX) {
151  VectorDouble &coords = getCoords();
152  double *coords_ptr = &*coords.data().begin();
153  size_t nb_integration_pts = getGaussPts().size2();
154  invJac.resize(4, nb_integration_pts, false);
156  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
157  double *ksi_ptr = &getGaussPts()(0, 0);
158  double *zeta_ptr = &getGaussPts()(1, 0);
159  for (size_t gg = 0; gg != nb_integration_pts;
160  ++gg, ++t_inv_jac, ++ksi_ptr, ++zeta_ptr) {
161  const double &ksi = *ksi_ptr;
162  const double &zeta = *zeta_ptr;
163  double j00 = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0x(zeta) +
164  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1x(zeta) +
165  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2x(zeta) +
166  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3x(zeta);
167  double j01 = coords_ptr[3 * 0 + 0] * diffN_MBQUAD0y(ksi) +
168  coords_ptr[3 * 1 + 0] * diffN_MBQUAD1y(ksi) +
169  coords_ptr[3 * 2 + 0] * diffN_MBQUAD2y(ksi) +
170  coords_ptr[3 * 3 + 0] * diffN_MBQUAD3y(ksi);
171  double j10 = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0x(zeta) +
172  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1x(zeta) +
173  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2x(zeta) +
174  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3x(zeta);
175  double j11 = coords_ptr[3 * 0 + 1] * diffN_MBQUAD0y(ksi) +
176  coords_ptr[3 * 1 + 1] * diffN_MBQUAD1y(ksi) +
177  coords_ptr[3 * 2 + 1] * diffN_MBQUAD2y(ksi) +
178  coords_ptr[3 * 3 + 1] * diffN_MBQUAD3y(ksi);
179  double det = j00 * j11 - j01 * j10;
180  t_inv_jac(0, 0) = j11 / det;
181  t_inv_jac(0, 1) = -j01 / det;
182  t_inv_jac(1, 0) = -j10 / det;
183  t_inv_jac(1, 1) = j00 / det;
184  }
185  }
187  };
188 
189  switch (getNumeredEntFiniteElementPtr()->getEntType()) {
190  case MBTRI:
191  CHKERR cal_inv_jac_on_tri();
192  break;
193  case MBQUAD:
194  CHKERR cal_inv_jac_on_quad();
195  break;
196  default:
197  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
198  "Operator not implemented for this entity type");
199  };
200 
201  doEntities[MBVERTEX] = true;
202  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
203 
205 }
206 
208 OpSetInvJacH1ForFace::doWork(int side, EntityType type,
211 
212  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI &&
213  getNumeredEntFiniteElementPtr()->getEntType() != MBQUAD)
214  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
215  "This operator can be used only with element which is triangle");
216 
217  auto apply_transform = [&](MatrixDouble &diff_n) {
219  size_t nb_functions = diff_n.size2() / 2;
220  if (nb_functions) {
221  size_t nb_gauss_pts = diff_n.size1();
222  diffNinvJac.resize(nb_gauss_pts, 2 * nb_functions, false);
223 
224  switch (type) {
225  case MBVERTEX:
226  case MBEDGE:
227  case MBTRI:
228  case MBQUAD: {
233  &diffNinvJac(0, 0), &diffNinvJac(0, 1));
235  &diff_n(0, 0), &diff_n(0, 1));
237  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
238  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
239  for (size_t dd = 0; dd != nb_functions; ++dd) {
240  t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
241  ++t_diff_n;
242  ++t_diff_n_ref;
243  }
244  }
245  diff_n.data().swap(diffNinvJac.data());
246  } break;
247  default:
248  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
249  }
250  }
252  };
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 
271 
273 }
274 
276 OpSetInvJacHcurlFace::doWork(int side, EntityType type,
279 
280  if (type != MBEDGE && type != MBTRI)
282 
283  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI)
284  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
285  "This operator can be used only with element which is triangle");
286 
290 
291  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
292 
293  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
294 
295  const unsigned int nb_base_functions = data.getDiffN(base).size2() / 6;
296  if (nb_base_functions) {
297  const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
298 
299  diffHcurlInvJac.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
300 
301  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
302  double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
304  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
305 
306  &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
307 
308  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1]);
309 
311  &invJac(0, 0), &invJac(1, 0), &invJac(2, 0), &invJac(3, 0));
312 
313  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
314  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
315  t_inv_diff_n(i, j) = t_diff_n(i, k) * t_inv_jac(k, j);
316  ++t_diff_n;
317  ++t_inv_diff_n;
318  }
319  }
320 
321  data.getDiffN(base).data().swap(diffHcurlInvJac.data());
322  }
323  }
324 
326 }
327 
329 OpMakeHdivFromHcurl::doWork(int side, EntityType type,
332 
333  if (type != MBEDGE && type != MBTRI)
335 
336  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI)
337  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
338  "This operator can be used only with element which is triangle");
339 
340  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
341 
342  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
343 
344  const size_t nb_base_functions = data.getN(base).size2() / 3;
345  if (nb_base_functions) {
346 
347  const size_t nb_gauss_pts = data.getN(base).size1();
348 
349  auto t_n = data.getFTensor1N<3>(base);
350  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
351  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
352  for (size_t bb = 0; bb != nb_base_functions; ++bb) {
353 
354  const double a = t_n(0);
355  t_n(0) = -t_n(1);
356  t_n(1) = a;
357 
358  for (auto n : {0, 1}) {
359  const double b = t_diff_n(0, n);
360  t_diff_n(0, n) = -t_diff_n(1, n);
361  t_diff_n(1, n) = b;
362  }
363 
364  ++t_n;
365  ++t_diff_n;
366  }
367  }
368  }
369  }
370 
372 }
373 
375  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
376 
378 
379  if (type != MBEDGE && type != MBTRI)
381 
385 
386  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
387 
388  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
389 
390  const size_t nb_base_functions = data.getN(base).size2() / 3;
391  if (nb_base_functions) {
392 
393  const size_t nb_gauss_pts = data.getN(base).size1();
394  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
395  if (data.getN(base).size2() > 0) {
396  auto t_n = data.getFTensor1N<3>(base);
397  double *t_transformed_n_ptr = &*piolaN.data().begin();
399  t_transformed_n_ptr, // HVEC0
400  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
402  &jAc(0, 0), &jAc(1, 0), &jAc(2, 0), &jAc(3, 0));
403  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
404  double det;
405  CHKERR determinantTensor2by2(t_jac, det);
406  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
407  t_transformed_n(i) = t_jac(i, k) * t_n(k) / det;
408  ++t_n;
409  ++t_transformed_n;
410  }
411  }
412  data.getN(base).data().swap(piolaN.data());
413  }
414 
415  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
416  if (data.getDiffN(base).size2() > 0) {
417  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
418  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
420  t_transformed_diff_n(t_transformed_diff_n_ptr,
421  &t_transformed_diff_n_ptr[HVEC0_1],
422 
423  &t_transformed_diff_n_ptr[HVEC1_0],
424  &t_transformed_diff_n_ptr[HVEC1_1],
425 
426  &t_transformed_diff_n_ptr[HVEC2_0],
427  &t_transformed_diff_n_ptr[HVEC2_1]);
429  &jAc(0, 0), &jAc(1, 0), &jAc(2, 0), &jAc(3, 0));
430  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
431  double det;
432  CHKERR determinantTensor2by2(t_jac, det);
433  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
434  t_transformed_diff_n(i, k) = t_jac(i, j) * t_diff_n(j, k) / det;
435  ++t_diff_n;
436  ++t_transformed_diff_n;
437  }
438  }
439  data.getDiffN(base).data().swap(piolaDiffN.data());
440  }
441  }
442  }
443 
445 }
446 
448  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
450 
451  if (type != MBEDGE)
453 
454  const auto &edge_direction = getDirection();
455 
457  FTensor::Tensor1<double, 3> t_m(-edge_direction[1], edge_direction[0],
458  edge_direction[2]);
459  const double l0 = t_m(i) * t_m(i);
460 
461  std::vector<double> l1;
462  {
463  int nb_gauss_pts = getTangetAtGaussPts().size1();
464  if (nb_gauss_pts) {
465  l1.resize(nb_gauss_pts);
466  const auto &edge_direction = getTangetAtGaussPts();
468  &edge_direction(0, 0), &edge_direction(0, 1), &edge_direction(0, 2));
469  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
470  l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
471  ++t_m_at_pts;
472  }
473  }
474  }
475 
476  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
477 
478  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
479  size_t nb_gauss_pts = data.getN(base).size1();
480  size_t nb_dofs = data.getN(base).size2() / 3;
481  if (nb_gauss_pts > 0 && nb_dofs > 0) {
482 
483  auto t_h_div = data.getFTensor1N<3>(base);
484 
485  size_t cc = 0;
486  const auto &edge_direction_at_gauss_pts = getTangetAtGaussPts();
487  if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
488 
490  &edge_direction_at_gauss_pts(0, 1),
491  &edge_direction_at_gauss_pts(0, 0),
492  &edge_direction_at_gauss_pts(0, 2));
493 
494  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
495  const double l0 = l1[gg];
496  for (int ll = 0; ll != nb_dofs; ++ll) {
497  const double val = t_h_div(0);
498  const double a = val / l0;
499  t_h_div(i) = t_m_at_pts(i) * a;
500  t_h_div(0) *= -1;
501  ++t_h_div;
502  ++cc;
503  }
504  ++t_m_at_pts;
505  }
506 
507  } else {
508 
509  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
510  for (int ll = 0; ll != nb_dofs; ll++) {
511  const double val = t_h_div(0);
512  const double a = val / l0;
513  t_h_div(i) = t_m(i) * a;
514  ++t_h_div;
515  ++cc;
516  }
517  }
518  }
519 
520  if (cc != nb_gauss_pts * nb_dofs)
521  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Data inconsistency");
522  }
523  }
524 
526 }
527 
529  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
531 
532  if (type == MBVERTEX) {
533 
534  VectorDouble &coords = getCoords();
535  double *coords_ptr = &*coords.data().begin();
536 
537  const int nb_gauss_pts = data.getN(NOBASE).size1();
538  auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
539 
543 
544  auto t_w = getFTensor0IntegrationWeight();
545  for (int gg = 0; gg != nb_gauss_pts; gg++) {
546 
547  FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
548  &coords_ptr[2], 3);
549  t_jac(i, j) = 0;
550  for (int bb = 0; bb != 6; bb++) {
551  t_jac(i, j) += t_coords(i) * t_diff_n(j);
552  ++t_diff_n;
553  ++t_coords;
554  }
555 
556  double det;
557  CHKERR determinantTensor3by3(t_jac, det);
558  t_w *= det / 2.;
559 
560  ++t_w;
561  }
562 
563  double &vol = getVolume();
564  auto t_w_scaled = getFTensor0IntegrationWeight();
565  for (int gg = 0; gg != nb_gauss_pts; gg++) {
566  t_w_scaled /= vol;
567  ++t_w_scaled;
568  }
569 
570  }
571 
572  doEntities[MBVERTEX] = true;
573  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
574 
576 }
577 
582 
583  if (type == MBVERTEX) {
584 
585  VectorDouble &coords = getCoords();
586  double *coords_ptr = &*coords.data().begin();
587 
588  const int nb_gauss_pts = data.getN(NOBASE).size1();
589  auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
590  invJac.resize(9, nb_gauss_pts, false);
591  invJac.clear();
592  auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
593 
597 
598  auto t_w = getFTensor0IntegrationWeight();
599  for (int gg = 0; gg != nb_gauss_pts; gg++) {
600 
601  FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
602  &coords_ptr[2], 3);
603  t_jac(i, j) = 0;
604  for (int bb = 0; bb != 6; bb++) {
605  t_jac(i, j) += t_coords(i) * t_diff_n(j);
606  ++t_diff_n;
607  ++t_coords;
608  }
609 
610  double det;
611  CHKERR determinantTensor3by3(t_jac, det);
612  CHKERR invertTensor3by3(t_jac, det, t_inv_jac);
613  ++t_inv_jac;
614 
615  }
616  }
617 
618  doEntities[MBVERTEX] = true;
619  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
620 
622 }
623 
628 
629  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
630 
631  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
632  if (data.getN(base).size2() == 0)
633  continue;
634 
635  const int nb_gauss_pts = data.getN(base).size1();
636  auto t_diff_n = data.getFTensor1DiffN<3>(base);
637  diffNinvJac.resize(data.getDiffN(base).size1(), data.getDiffN(base).size2(),
638  false);
639 
642 
644  &diffNinvJac(0, 0), &diffNinvJac(0, 1), &diffNinvJac(0, 2));
645  auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
646 
647  const int nb_dofs = data.getN(base).size2();
648  for (int gg = 0; gg != nb_gauss_pts; gg++) {
649  for (int bb = 0; bb != nb_dofs; bb++) {
650  t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
651  ++t_inv_diff_n;
652  ++t_diff_n;
653  }
654  ++t_inv_jac;
655  }
656 
657  data.getDiffN(base).data().swap(diffNinvJac.data());
658  }
659 
661 }
662 
666 
668 
669  if (type == MBVERTEX) {
670 
671  VectorDouble &coords = getCoords();
672  double *coords_ptr = &*coords.data().begin();
673  double diff_n[6];
674  CHKERR ShapeDiffMBTRI(diff_n);
675  double j00_f3, j01_f3, j10_f3, j11_f3;
676  for (int gg = 0; gg < 1; gg++) {
677  // this is triangle, derivative of nodal shape functions is constant.
678  // So only need to do one node.
679  j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
680  j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
681  j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
682  j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
683  }
684  double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
685  invJacF3.resize(2, 2, false);
686  invJacF3(0, 0) = j11_f3 / det_f3;
687  invJacF3(0, 1) = -j01_f3 / det_f3;
688  invJacF3(1, 0) = -j10_f3 / det_f3;
689  invJacF3(1, 1) = j00_f3 / det_f3;
690  }
691 
692  doEntities[MBVERTEX] = true;
693  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
694 
696 }
697 
702 
703  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
704 
705  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
706 
707  unsigned int nb_dofs = data.getN(base).size2();
708  if (nb_dofs == 0)
710  unsigned int nb_gauss_pts = data.getN(base).size1();
711  diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs, false);
712 
713  if (type != MBVERTEX) {
714  if (nb_dofs != data.getDiffN(base).size2() / 2) {
715  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
716  "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
717  "%u/2 )",
718  nb_dofs, data.getDiffN(base).size2());
719  }
720  }
721 
722  switch (type) {
723  case MBVERTEX:
724  case MBEDGE:
725  case MBTRI: {
726  for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
727  for (unsigned int dd = 0; dd < nb_dofs; dd++) {
729  &*invJacF3.data().begin(), 2,
730  &data.getDiffN(base)(gg, 2 * dd), 1, 0,
731  &diffNinvJac(gg, 2 * dd), 1);
732  }
733  }
734  data.getDiffN(base).data().swap(diffNinvJac.data());
735  } break;
736  default:
737  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
738  }
739  }
740 
742 }
743 
744 } // 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:482
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:507
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:514
#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
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:495
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:463
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:150
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
auto getFTensor0IntegrationWeight()
Get integration weights.
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:602
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
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
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
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
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:72