v0.10.0
DataOperators.cpp
Go to the documentation of this file.
1 /** file DataOperators.cpp
2 
3  \brief implementation of Data Operators 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 #ifdef __cplusplus
22 extern "C" {
23 #endif
24 #include <cblas.h>
25 #include <lapack_wrap.h>
26 #include <gm_rule.h>
27 #ifdef __cplusplus
28 }
29 #endif
30 
31 namespace MoFEM {
32 
34  :
35 
36  sYmm(symm),
37 
38  doEntities{true, true, true, true, true, true,
39  true, true, true, true, true, true},
40 
41  doVertices(doEntities[MBVERTEX]), doEdges(doEntities[MBEDGE]),
42  doQuads(doEntities[MBQUAD]), doTris(doEntities[MBTRI]),
43  doTets(doEntities[MBTET]), doPrisms(doEntities[MBPRISM]) {
44 
45  /// This not yet implemented, switch off.
46  doEntities[MBPOLYGON] = false;
47  doEntities[MBPYRAMID] = false;
48  doEntities[MBKNIFE] = false;
49  doEntities[MBHEX] = false;
50  doEntities[MBPOLYHEDRON] = false;
51 }
52 
53 template <bool Symm>
55  DataForcesAndSourcesCore &col_data) {
57 
58  auto do_col_entity =
59  [&](boost::ptr_vector<DataForcesAndSourcesCore::EntData> &row_ent_data,
60  const int ss, const EntityType row_type, const EntityType low_type,
61  const EntityType hi_type) {
63  for (EntityType col_type = low_type; col_type != hi_type; ++col_type) {
64  auto &col_ent_data = col_data.dataOnEntities[col_type];
65  for (size_t SS = 0; SS != col_ent_data.size(); SS++) {
66  if (col_ent_data[SS].getFieldData().size())
67  CHKERR doWork(ss, SS, row_type, col_type, row_ent_data[ss],
68  col_ent_data[SS]);
69  }
70  }
72  };
73 
74  auto do_row_entity = [&](const EntityType type) {
76  auto &row_ent_data = row_data.dataOnEntities[type];
77  for (size_t ss = 0; ss != row_ent_data.size(); ++ss) {
78  size_t SS = 0;
79  if (Symm)
80  SS = ss;
81  for (; SS < col_data.dataOnEntities[type].size(); SS++) {
82  CHKERR doWork(ss, SS, type, type, row_ent_data[ss],
83  col_data.dataOnEntities[type][SS]);
84  }
85  if (!Symm)
86  CHKERR do_col_entity(row_ent_data, ss, type, MBVERTEX, type);
87  CHKERR do_col_entity(row_ent_data, ss, type,
88  static_cast<EntityType>(type + 1), MBMAXTYPE);
89  }
91  };
92 
93  for (EntityType row_type = MBVERTEX; row_type != MBENTITYSET; ++row_type) {
94  if (doEntities[row_type]) {
95  CHKERR do_row_entity(row_type);
96  }
97  }
98 
99  if (doEntities[MBENTITYSET]) {
100  for (unsigned int mm = 0; mm != row_data.dataOnEntities[MBENTITYSET].size();
101  ++mm) {
102  if (!row_data.dataOnEntities[MBENTITYSET][mm].getFieldData().empty()) {
103  CHKERR do_row_entity(MBENTITYSET);
104  }
105  }
106  }
107 
109 }
110 
112  DataForcesAndSourcesCore &col_data) {
113  if (getSymm())
114  return opLhs<true>(row_data, col_data);
115  else
116  return opLhs<false>(row_data, col_data);
117 }
118 
119 template <bool ErrorIfNoBase>
122  const std::array<bool, MBMAXTYPE> &do_entities) {
124 
125  auto do_entity = [&](auto type) {
127 
128  auto &ent_data = data.dataOnEntities[type];
129  const size_t size = ent_data.size();
130  for (size_t ss = 0; ss != size; ++ss) {
131 
132  auto &side_data = ent_data[ss];
133 
134  if (ErrorIfNoBase) {
135  if (side_data.getFieldData().size() &&
136  (side_data.getBase() == NOBASE ||
137  side_data.getBase() == LASTBASE)) {
138  for (VectorDofs::iterator it = side_data.getFieldDofs().begin();
139  it != side_data.getFieldDofs().end(); it++)
140  if ((*it) && (*it)->getActive())
141  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "No base on");
142  }
143  }
144 
145  CHKERR doWork(ss, type, side_data);
146  }
147 
149  };
150 
151  for (EntityType row_type = MBVERTEX; row_type != MBENTITYSET; ++row_type) {
152  if (do_entities[row_type]) {
153  CHKERR do_entity(row_type);
154  }
155  }
156 
157  if (do_entities[MBENTITYSET]) {
158  // This is odd behaviour, diffrent than for other entities. Should be
159  // changed that behaviour is consistent,
160  for (unsigned int mm = 0; mm != data.dataOnEntities[MBENTITYSET].size();
161  ++mm) {
162  if (!data.dataOnEntities[MBENTITYSET][mm].getFieldData().empty()) {
163  CHKERR doWork(mm, MBENTITYSET, data.dataOnEntities[MBENTITYSET][mm]);
164  }
165  }
166  }
167 
169 }
170 
172  const bool error_if_no_base) {
173  if (error_if_no_base)
174  return opRhs<true>(data, doEntities);
175  else
176  return opRhs<false>(data, doEntities);
177 }
178 
179 template <>
180 MoFEMErrorCode invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
181  MatrixDouble &jac_data, VectorDouble &det_data,
182  MatrixDouble &inv_jac_data) {
184  auto A = getFTensor2FromMat<3, 3>(jac_data);
185  int nb_gauss_pts = jac_data.size2();
186  det_data.resize(nb_gauss_pts, false);
187  inv_jac_data.resize(3, nb_gauss_pts, false);
188  auto det = getFTensor0FromVec(det_data);
189  auto I = getFTensor2FromMat<3, 3>(inv_jac_data);
190  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
192  CHKERR invertTensor3by3(A, det, I);
193  ++A;
194  ++det;
195  ++I;
196  }
198 }
199 
203 
204  auto transform_base = [&](MatrixDouble &diff_n,
205  const bool diff_at_gauss_ptr) {
207 
208  if (!diff_n.size1())
210  if (!diff_n.size2())
212 
213  const int nb_base_functions =
214  (diff_at_gauss_ptr || type != MBVERTEX) ? diff_n.size2() / 3 : 4;
215  const int nb_gauss_pts =
216  (diff_at_gauss_ptr || type != MBVERTEX) ? diff_n.size1() : 1;
217  diffNinvJac.resize(diff_n.size1(), diff_n.size2(), false);
218 
219  double *t_diff_n_ptr = &*diff_n.data().begin();
221  t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
222  double *t_inv_n_ptr = &*diffNinvJac.data().begin();
224  t_inv_n_ptr, &t_inv_n_ptr[1], &t_inv_n_ptr[2]);
225 
226  switch (type) {
227 
228  case MBVERTEX:
229  case MBEDGE:
230  case MBTRI:
231  case MBTET: {
232  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
233  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
234  t_inv_diff_n(i) = t_diff_n(j) * tInvJac(j, i);
235  ++t_diff_n;
236  ++t_inv_diff_n;
237  }
238  }
239 
240  } break;
241  default:
242  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
243  }
244 
245  diff_n.data().swap(diffNinvJac.data());
246 
248  };
249 
250  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
251  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
252  CHKERR transform_base(data.getDiffN(base), false);
253  }
254 
255  switch (type) {
256  case MBVERTEX:
257  for (auto &m : data.getBBDiffNMap())
258  CHKERR transform_base(*(m.second), true);
259  break;
260  default:
261  for (auto &ptr : data.getBBDiffNByOrderArray())
262  if (ptr)
263  CHKERR transform_base(*ptr, true);
264  }
265 
267 }
268 
273 
274  if (type != MBEDGE && type != MBTRI && type != MBTET)
276 
277  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
278 
279  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
280 
281  const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
282  const unsigned int nb_base_functions = data.getDiffN(base).size2() / 9;
283  if (!nb_base_functions)
284  continue;
285 
286  diffHdivInvJac.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
287 
288  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
289  double *t_inv_diff_n_ptr = &*diffHdivInvJac.data().begin();
291  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
292  &t_inv_diff_n_ptr[HVEC0_2],
293 
294  &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
295  &t_inv_diff_n_ptr[HVEC1_2],
296 
297  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
298  &t_inv_diff_n_ptr[HVEC2_2]);
299 
300  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
301  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
302  t_inv_diff_n(k, i) = t_diff_n(k, j) * tInvJac(j, i);
303  ++t_diff_n;
304  ++t_inv_diff_n;
305  }
306  }
307 
308  data.getDiffN(base).data().swap(diffHdivInvJac.data());
309  }
310 
312 }
313 
315  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
317 
318  if (type != MBTRI && type != MBTET)
320 
321  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
322 
323  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
324 
325  const unsigned int nb_base_functions = data.getN(base).size2() / 3;
326  if (!nb_base_functions)
327  continue;
328 
329  const double c = 1. / 6.;
330  const unsigned int nb_gauss_pts = data.getN(base).size1();
331 
332  double const a = c / vOlume;
333 
334  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
335  if (data.getN(base).size2() > 0) {
336  auto t_n = data.getFTensor1N<3>(base);
337  double *t_transformed_n_ptr = &*piolaN.data().begin();
339  t_transformed_n_ptr, // HVEC0
340  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
341  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
342  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
343  t_transformed_n(i) = a * tJac(i, k) * t_n(k);
344  ++t_n;
345  ++t_transformed_n;
346  }
347  }
348  data.getN(base).data().swap(piolaN.data());
349  }
350 
351  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
352  if (data.getDiffN(base).size2() > 0) {
353  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
354  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
356  t_transformed_diff_n(t_transformed_diff_n_ptr,
357  &t_transformed_diff_n_ptr[HVEC0_1],
358  &t_transformed_diff_n_ptr[HVEC0_2],
359  &t_transformed_diff_n_ptr[HVEC1_0],
360  &t_transformed_diff_n_ptr[HVEC1_1],
361  &t_transformed_diff_n_ptr[HVEC1_2],
362  &t_transformed_diff_n_ptr[HVEC2_0],
363  &t_transformed_diff_n_ptr[HVEC2_1],
364  &t_transformed_diff_n_ptr[HVEC2_2]);
365  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
366  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
367  t_transformed_diff_n(i, k) = a * tJac(i, j) * t_diff_n(j, k);
368  ++t_diff_n;
369  ++t_transformed_diff_n;
370  }
371  }
372  data.getDiffN(base).data().swap(piolaDiffN.data());
373  }
374  }
375 
377 }
378 
383 
384  if (type != MBEDGE && type != MBTRI && type != MBTET)
386 
387  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
388 
389  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
390 
391  const unsigned int nb_base_functions = data.getN(base).size2() / 3;
392  if (!nb_base_functions)
393  continue;
394 
395  const unsigned int nb_gauss_pts = data.getN(base).size1();
396  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
397  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
398 
399  auto t_n = data.getFTensor1N<3>(base);
400  double *t_transformed_n_ptr = &*piolaN.data().begin();
402  t_transformed_n_ptr, &t_transformed_n_ptr[HVEC1],
403  &t_transformed_n_ptr[HVEC2]);
404  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
405  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
406  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
407  t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
408  &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
409  &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
410  &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
411  &t_transformed_diff_n_ptr[HVEC2_2]);
412 
413  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
414  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
415  t_transformed_n(i) = tInvJac(k, i) * t_n(k);
416  t_transformed_diff_n(i, k) = tInvJac(j, i) * t_diff_n(j, k);
417  ++t_n;
418  ++t_transformed_n;
419  ++t_diff_n;
420  ++t_transformed_diff_n;
421  }
422  }
423  data.getN(base).data().swap(piolaN.data());
424  data.getDiffN(base).data().swap(piolaDiffN.data());
425  }
426 
427  // data.getBase() = base;
428 
430 }
431 
433 OpSetHoInvJacH1::doWork(int side, EntityType type,
436 
437  if (invHoJac.size2() != 9)
438  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
439  "It looks that ho inverse of Jacobian is not calculated %d != 9",
440  invHoJac.size2());
441 
442  auto transform_base = [&](MatrixDouble &diff_n) {
444 
445  unsigned int nb_gauss_pts = diff_n.size1();
446  if (nb_gauss_pts == 0)
448 
449  if (invHoJac.size1() == nb_gauss_pts) {
450 
451  unsigned int nb_base_functions = diff_n.size2() / 3;
452  if (nb_base_functions == 0)
454 
455  double *t_inv_jac_ptr = &*invHoJac.data().begin();
457  t_inv_jac_ptr, &t_inv_jac_ptr[1], &t_inv_jac_ptr[2],
458  &t_inv_jac_ptr[3], &t_inv_jac_ptr[4], &t_inv_jac_ptr[5],
459  &t_inv_jac_ptr[6], &t_inv_jac_ptr[7], &t_inv_jac_ptr[8], 9);
460 
461  diffNinvJac.resize(nb_gauss_pts, 3 * nb_base_functions, false);
462 
463  double *t_diff_n_ptr = &*diff_n.data().begin();
465  t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
466  double *t_inv_n_ptr = &*diffNinvJac.data().begin();
467  FTensor::Tensor1<double *, 3> t_inv_diff_n(t_inv_n_ptr, &t_inv_n_ptr[1],
468  &t_inv_n_ptr[2], 3);
469 
470  switch (type) {
471  case MBVERTEX:
472  case MBEDGE:
473  case MBTRI:
474  case MBTET: {
475  for (unsigned int gg = 0; gg < nb_gauss_pts; ++gg) {
476  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
477  t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
478  ++t_diff_n;
479  ++t_inv_diff_n;
480  }
481  ++t_inv_jac;
482  }
483  } break;
484  default:
485  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
486  }
487 
488  diff_n.data().swap(diffNinvJac.data());
489  }
491  };
492 
493  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
494  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
495  CHKERR transform_base(data.getDiffN(base));
496  }
497 
498  switch (type) {
499  case MBVERTEX:
500  for (auto &m : data.getBBDiffNMap())
501  CHKERR transform_base(*(m.second));
502  break;
503  default:
504  for (auto &ptr : data.getBBDiffNByOrderArray())
505  if (ptr)
506  CHKERR transform_base(*ptr);
507  }
508 
510 }
511 
516 
517  if (type != MBEDGE && type != MBTRI && type != MBTET)
519 
520  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
521 
522  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
523 
524  diffHdivInvJac.resize(data.getDiffN(base).size1(),
525  data.getDiffN(base).size2(), false);
526 
527  unsigned int nb_gauss_pts = data.getDiffN(base).size1();
528  unsigned int nb_base_functions = data.getDiffN(base).size2() / 9;
529  if (nb_base_functions == 0)
530  continue;
531 
532  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
533  double *t_inv_diff_n_ptr = &*diffHdivInvJac.data().begin();
535  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
536  &t_inv_diff_n_ptr[HVEC0_2], &t_inv_diff_n_ptr[HVEC1_0],
537  &t_inv_diff_n_ptr[HVEC1_1], &t_inv_diff_n_ptr[HVEC1_2],
538  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
539  &t_inv_diff_n_ptr[HVEC2_2]);
540  double *t_inv_jac_ptr = &*invHoJac.data().begin();
542  t_inv_jac_ptr, &t_inv_jac_ptr[1], &t_inv_jac_ptr[2], &t_inv_jac_ptr[3],
543  &t_inv_jac_ptr[4], &t_inv_jac_ptr[5], &t_inv_jac_ptr[6],
544  &t_inv_jac_ptr[7], &t_inv_jac_ptr[8]);
545 
546  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
547  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
548  t_inv_diff_n(i, j) = t_inv_jac(k, j) * t_diff_n(i, k);
549  ++t_diff_n;
550  ++t_inv_diff_n;
551  }
552  ++t_inv_jac;
553  }
554 
555  data.getDiffN(base).data().swap(diffHdivInvJac.data());
556  }
557 
559 }
560 
562  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
564 
565  if (type != MBTRI && type != MBTET)
567 
568  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
569 
570  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
571 
572  unsigned int nb_gauss_pts = data.getN(base).size1();
573  unsigned int nb_base_functions = data.getN(base).size2() / 3;
574  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
575  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
576 
577  auto t_n = data.getFTensor1N<3>(base);
578  double *t_transformed_n_ptr = &*piolaN.data().begin();
580  t_transformed_n_ptr, // HVEC0
581  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
582  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
583  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
584  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
585  t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
586  &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
587  &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
588  &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
589  &t_transformed_diff_n_ptr[HVEC2_2]);
590 
591  FTensor::Tensor0<double *> t_det(&*detHoJac.data().begin());
592  double *t_jac_ptr = &*hoJac.data().begin();
594  t_jac_ptr, &t_jac_ptr[1], &t_jac_ptr[2], &t_jac_ptr[3], &t_jac_ptr[4],
595  &t_jac_ptr[5], &t_jac_ptr[6], &t_jac_ptr[7], &t_jac_ptr[8]);
596 
597  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
598  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
599  const double a = 1. / t_det;
600  t_transformed_n(i) = a * t_jac(i, k) * t_n(k);
601  t_transformed_diff_n(i, k) = a * t_jac(i, j) * t_diff_n(j, k);
602  ++t_n;
603  ++t_transformed_n;
604  ++t_diff_n;
605  ++t_transformed_diff_n;
606  }
607  ++t_det;
608  ++t_jac;
609  }
610 
611  data.getN(base).data().swap(piolaN.data());
612  data.getDiffN(base).data().swap(piolaDiffN.data());
613  }
614 
616 }
617 
619  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
621 
622  if (type != MBEDGE && type != MBTRI && type != MBTET)
624 
625  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
626 
627  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
628 
629  unsigned int nb_gauss_pts = data.getN(base).size1();
630  unsigned int nb_base_functions = data.getN(base).size2() / 3;
631  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
632  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
633 
634  auto t_n = data.getFTensor1N<3>(base);
635  double *t_transformed_n_ptr = &*piolaN.data().begin();
637  t_transformed_n_ptr, // HVEC0
638  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
639  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
640  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
641  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
642  t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
643  &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
644  &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
645  &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
646  &t_transformed_diff_n_ptr[HVEC2_2]);
647 
648  double *t_inv_jac_ptr = &*hoInvJac.data().begin();
650  t_inv_jac_ptr, &t_inv_jac_ptr[1], &t_inv_jac_ptr[2], &t_inv_jac_ptr[3],
651  &t_inv_jac_ptr[4], &t_inv_jac_ptr[5], &t_inv_jac_ptr[6],
652  &t_inv_jac_ptr[7], &t_inv_jac_ptr[8], 9);
653 
654  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
655  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
656  t_transformed_n(i) = t_inv_jac(k, i) * t_n(k);
657  t_transformed_diff_n(i, k) = t_inv_jac(j, i) * t_diff_n(j, k);
658  ++t_n;
659  ++t_transformed_n;
660  ++t_diff_n;
661  ++t_transformed_diff_n;
662  }
663  ++t_inv_jac;
664  }
665 
666  data.getN(base).data().swap(piolaN.data());
667  data.getDiffN(base).data().swap(piolaDiffN.data());
668  }
669 
671 }
672 
677 
678  unsigned int nb_dofs = data.getFieldData().size();
679  if (nb_dofs == 0)
681 
682  int nb_gauss_pts = data.getN().size1();
683  cOords_at_GaussPt.resize(nb_gauss_pts, 3, false);
684  tAngent1_at_GaussPt.resize(nb_gauss_pts, 3, false);
685  tAngent2_at_GaussPt.resize(nb_gauss_pts, 3, false);
686 
687  auto get_ftensor1 = [](MatrixDouble &m) {
689  &m(0, 0), &m(0, 1), &m(0, 2));
690  };
691  auto t_coords = get_ftensor1(cOords_at_GaussPt);
692  auto t_t1 = get_ftensor1(tAngent1_at_GaussPt);
693  auto t_t2 = get_ftensor1(tAngent2_at_GaussPt);
697 
698  switch (type) {
699  case MBVERTEX: {
700  cOords_at_GaussPt.clear();
701  tAngent1_at_GaussPt.clear();
702  tAngent2_at_GaussPt.clear();
703  }
704  case MBEDGE:
705  case MBTRI:
706  case MBQUAD: {
707  if (2 * data.getN().size2() != data.getDiffN().size2()) {
708  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
709  }
710  if (nb_dofs % 3 != 0) {
711  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
712  }
713  if (nb_dofs > 3 * data.getN().size2()) {
714  unsigned int nn = 0;
715  for (; nn != nb_dofs; nn++) {
716  if (!data.getFieldDofs()[nn]->getActive())
717  break;
718  }
719  if (nn > 3 * data.getN().size2()) {
720  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
721  "data inconsistency for base %s",
723  } else {
724  nb_dofs = nn;
725  if (!nb_dofs)
727  }
728  }
729  const int nb_base_functions = data.getN().size2();
730  auto t_base = data.getFTensor0N();
731  auto t_diff_base = data.getFTensor1DiffN<2>();
732  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
733  auto t_data = data.getFTensor1FieldData<3>();
734  int bb = 0;
735  for (; bb != nb_dofs / 3; ++bb) {
736  t_coords(i) += t_base * t_data(i);
737  t_t1(i) += t_data(i) * t_diff_base(N0);
738  t_t2(i) += t_data(i) * t_diff_base(N1);
739  ++t_data;
740  ++t_base;
741  ++t_diff_base;
742  }
743  for (; bb != nb_base_functions; ++bb) {
744  ++t_base;
745  ++t_diff_base;
746  }
747  ++t_coords;
748  ++t_t1;
749  ++t_t2;
750  }
751  } break;
752  default:
753  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
754  }
755 
757 }
758 
761 
762  nOrmals_at_GaussPt.resize(tAngent1_at_GaussPt.size1(), 3, false);
763 
764  auto get_ftensor1 = [](MatrixDouble &m) {
766  &m(0, 0), &m(0, 1), &m(0, 2));
767  };
768  auto t_normal = get_ftensor1(nOrmals_at_GaussPt);
769  auto t_t1 = get_ftensor1(tAngent1_at_GaussPt);
770  auto t_t2 = get_ftensor1(tAngent2_at_GaussPt);
771 
775 
776  for (unsigned int gg = 0; gg != tAngent1_at_GaussPt.size1(); ++gg) {
777  t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
778  ++t_normal;
779  ++t_t1;
780  ++t_t2;
781  }
782 
784 }
785 
790 
791  if (data.getFieldData().size() == 0)
793  const int valid_edges3[] = {1, 1, 1, 0, 0, 0, 0, 0, 0};
794  const int valid_faces3[] = {0, 0, 0, 1, 0, 0, 0, 0, 0};
795  const int valid_edges4[] = {0, 0, 0, 0, 0, 0, 1, 1, 1};
796  const int valid_faces4[] = {0, 0, 0, 0, 1, 0, 0, 0, 0};
797 
798  if (type == MBEDGE) {
799  if (!valid_edges3[side] || valid_edges4[side])
801  } else if (type == MBTRI) {
802  if (!valid_faces3[side] || valid_faces4[side])
804  }
805 
806  switch (type) {
807  case MBVERTEX: {
808  for (unsigned int gg = 0; gg < data.getN().size1(); ++gg) {
809  for (int dd = 0; dd < 3; dd++) {
810  cOords_at_GaussPtF3(gg, dd) =
811  cblas_ddot(3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
813  3, &data.getDiffN()(gg, 0), 2, &data.getFieldData()[dd], 3);
815  3, &data.getDiffN()(gg, 1), 2, &data.getFieldData()[dd], 3);
817  3, &data.getN(gg)[0], 1, &data.getFieldData()[9 + dd], 3);
819  3, &data.getDiffN()(gg, 6 + 0), 2, &data.getFieldData()[9 + dd], 3);
821  3, &data.getDiffN()(gg, 6 + 1), 2, &data.getFieldData()[9 + dd], 3);
822  }
823  }
824  } break;
825  case MBEDGE:
826  case MBTRI: {
827  if (2 * data.getN().size2() != data.getDiffN().size2()) {
828  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
829  }
830  unsigned int nb_dofs = data.getFieldData().size();
831  if (nb_dofs % 3 != 0) {
832  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
833  }
834  if (nb_dofs > 3 * data.getN().size2()) {
835  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
836  "data inconsistency, side %d type %d", side, type);
837  }
838  for (unsigned int gg = 0; gg < data.getN().size1(); ++gg) {
839  for (int dd = 0; dd < 3; dd++) {
840  if ((type == MBTRI && valid_faces3[side]) ||
841  (type == MBEDGE && valid_edges3[side])) {
843  nb_dofs / 3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
844  tAngent1_at_GaussPtF3(gg, dd) +=
845  cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 0), 2,
846  &data.getFieldData()[dd], 3);
847  tAngent2_at_GaussPtF3(gg, dd) +=
848  cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 1), 2,
849  &data.getFieldData()[dd], 3);
850  } else if ((type == MBTRI && valid_faces4[side]) ||
851  (type == MBEDGE && valid_edges4[side])) {
853  nb_dofs / 3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
854  tAngent1_at_GaussPtF4(gg, dd) +=
855  cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 0), 2,
856  &data.getFieldData()[dd], 3);
857  tAngent2_at_GaussPtF4(gg, dd) +=
858  cblas_ddot(nb_dofs / 3, &data.getDiffN()(gg, 1), 2,
859  &data.getFieldData()[dd], 3);
860  }
861  }
862  }
863  } break;
864  default:
865  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
866  }
867 
869 }
870 
873 
874  sPin.resize(3, 3);
875  sPin.clear();
876  nOrmals_at_GaussPtF3.resize(tAngent1_at_GaussPtF3.size1(), 3, false);
877  for (unsigned int gg = 0; gg < tAngent1_at_GaussPtF3.size1(); ++gg) {
878  ierr = Spin(&*sPin.data().begin(), &tAngent1_at_GaussPtF3(gg, 0));
879  CHKERRG(ierr);
880  cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., &*sPin.data().begin(), 3,
881  &tAngent2_at_GaussPtF3(gg, 0), 1, 0.,
882  &nOrmals_at_GaussPtF3(gg, 0), 1);
883  }
884  sPin.clear();
885  nOrmals_at_GaussPtF4.resize(tAngent1_at_GaussPtF4.size1(), 3, false);
886  for (unsigned int gg = 0; gg < tAngent1_at_GaussPtF4.size1(); ++gg) {
887  ierr = Spin(&*sPin.data().begin(), &tAngent1_at_GaussPtF4(gg, 0));
888  CHKERRG(ierr);
889  cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., &*sPin.data().begin(), 3,
890  &tAngent2_at_GaussPtF4(gg, 0), 1, 0.,
891  &nOrmals_at_GaussPtF4(gg, 0), 1);
892  }
894 }
895 
897  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
900 
901  if (type != MBTRI)
903 
904  if (normalRawPtr == nullptr && normalsAtGaussPtsRawPtr == nullptr)
905  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
906  "Pointer to normal/normals not set");
907 
908  bool normal_is_at_gauss_pts = (normalsAtGaussPtsRawPtr != nullptr);
909  if (normal_is_at_gauss_pts)
910  normal_is_at_gauss_pts = (normalsAtGaussPtsRawPtr->size1() != 0);
911 
912  auto apply_transform_linear_geometry = [&](auto base, auto nb_gauss_pts,
913  auto nb_base_functions) {
915  const auto &normal = *normalRawPtr;
916  auto t_normal = FTensor::Tensor1<double, 3>{normal[normalShift + 0],
917  normal[normalShift + 1],
918  normal[normalShift + 2]};
919  const auto l02 = t_normal(i) * t_normal(i);
920  auto t_base = data.getFTensor1N<3>(base);
921  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
922  for (int bb = 0; bb != nb_base_functions; ++bb) {
923  const auto v = t_base(0);
924  t_base(i) = (v / l02) * t_normal(i);
925  ++t_base;
926  }
927  }
929  };
930 
931  auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
932  auto nb_base_functions) {
934  const MatrixDouble &normals_at_pts = *normalsAtGaussPtsRawPtr;
936  &normals_at_pts(0, 0), &normals_at_pts(0, 1), &normals_at_pts(0, 2));
937  auto t_base = data.getFTensor1N<3>(base);
938  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
939  for (int bb = 0; bb != nb_base_functions; ++bb) {
940  const auto v = t_base(0);
941  const auto l2 = t_normal(i) * t_normal(i);
942  t_base(i) = (v / l2) * t_normal(i);
943  ++t_base;
944  }
945  ++t_normal;
946  }
948  };
949 
950  if (normal_is_at_gauss_pts) {
951  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
952 
953  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
954  const auto &base_functions = data.getN(base);
955  const auto nb_gauss_pts = base_functions.size1();
956 
957  if (nb_gauss_pts) {
958 
959  if (normalsAtGaussPtsRawPtr->size1() != nb_gauss_pts)
960  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
961  "normalsAtGaussPtsRawPtr has inconsistent number of "
962  "integration "
963  "points");
964 
965  const auto nb_base_functions = base_functions.size2() / 3;
966  CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
967  nb_base_functions);
968  }
969  }
970  } else {
971  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
972 
973  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
974  const auto &base_functions = data.getN(base);
975  const auto nb_gauss_pts = base_functions.size1();
976 
977  if (nb_gauss_pts) {
978  const auto nb_base_functions = base_functions.size2() / 3;
979  CHKERR apply_transform_linear_geometry(base, nb_gauss_pts,
980  nb_base_functions);
981  }
982  }
983  }
984 
986 }
987 
989  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
991 
992  if (type != MBEDGE && type != MBTRI)
994 
998 
1000  &tAngent0[0], &tAngent1[0], &nOrmal[0],
1001 
1002  &tAngent0[1], &tAngent1[1], &nOrmal[1],
1003 
1004  &tAngent0[2], &tAngent1[2], &nOrmal[2],
1005 
1006  3);
1007  double det;
1009  CHKERR determinantTensor3by3(t_m, det);
1010  CHKERR invertTensor3by3(t_m, det, t_inv_m);
1011 
1012  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
1013 
1014  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
1015 
1016  auto &baseN = data.getN(base);
1017  auto &diffBaseN = data.getDiffN(base);
1018 
1019  int nb_dofs = baseN.size2() / 3;
1020  int nb_gauss_pts = baseN.size1();
1021 
1022  MatrixDouble piola_n(baseN.size1(), baseN.size2());
1023  MatrixDouble diff_piola_n(diffBaseN.size1(), diffBaseN.size2());
1024 
1025  if (nb_dofs > 0 && nb_gauss_pts > 0) {
1026 
1028  &baseN(0, HVEC0), &baseN(0, HVEC1), &baseN(0, HVEC2));
1030  &diffBaseN(0, HVEC0_0), &diffBaseN(0, HVEC0_1),
1031  &diffBaseN(0, HVEC1_0), &diffBaseN(0, HVEC1_1),
1032  &diffBaseN(0, HVEC2_0), &diffBaseN(0, HVEC2_1));
1033  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
1034  &piola_n(0, HVEC0), &piola_n(0, HVEC1), &piola_n(0, HVEC2));
1036  t_transformed_diff_h_curl(
1037  &diff_piola_n(0, HVEC0_0), &diff_piola_n(0, HVEC0_1),
1038  &diff_piola_n(0, HVEC1_0), &diff_piola_n(0, HVEC1_1),
1039  &diff_piola_n(0, HVEC2_0), &diff_piola_n(0, HVEC2_1));
1040 
1041  int cc = 0;
1042  if (normalsAtGaussPts.size1() == (unsigned int)nb_gauss_pts) {
1043  // HO geometry is set, so jacobian is different at each gauss point
1045  &tangent0AtGaussPt(0, 0), &tangent1AtGaussPt(0, 0),
1046  &normalsAtGaussPts(0, 0), &tangent0AtGaussPt(0, 1),
1047  &tangent1AtGaussPt(0, 1), &normalsAtGaussPts(0, 1),
1048  &tangent0AtGaussPt(0, 2), &tangent1AtGaussPt(0, 2),
1049  &normalsAtGaussPts(0, 2), 3);
1050  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1051  CHKERR determinantTensor3by3(t_m_at_pts, det);
1052  CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
1053  for (int ll = 0; ll != nb_dofs; ll++) {
1054  t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
1055  t_transformed_diff_h_curl(i, k) =
1056  t_inv_m(j, i) * t_diff_h_curl(j, k);
1057  ++t_h_curl;
1058  ++t_transformed_h_curl;
1059  ++t_diff_h_curl;
1060  ++t_transformed_diff_h_curl;
1061  ++cc;
1062  }
1063  ++t_m_at_pts;
1064  }
1065  } else {
1066  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1067  for (int ll = 0; ll != nb_dofs; ll++) {
1068  t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
1069  t_transformed_diff_h_curl(i, k) =
1070  t_inv_m(j, i) * t_diff_h_curl(j, k);
1071  ++t_h_curl;
1072  ++t_transformed_h_curl;
1073  ++t_diff_h_curl;
1074  ++t_transformed_diff_h_curl;
1075  ++cc;
1076  }
1077  }
1078  }
1079  if (cc != nb_gauss_pts * nb_dofs)
1080  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Data inconsistency");
1081 
1082  baseN.data().swap(piola_n.data());
1083  diffBaseN.data().swap(diff_piola_n.data());
1084  }
1085  }
1086 
1088 }
1089 
1091 OpGetHoTangentOnEdge::doWork(int side, EntityType type,
1094 
1095  int nb_dofs = data.getFieldData().size();
1096  if (nb_dofs == 0)
1098 
1099  int nb_gauss_pts = data.getN().size1();
1100  tAngent.resize(nb_gauss_pts, 3, false);
1101 
1102  int nb_approx_fun = data.getN().size2();
1103  double *diff = &*data.getDiffN().data().begin();
1104  double *dofs[] = {&data.getFieldData()[0], &data.getFieldData()[1],
1105  &data.getFieldData()[2]};
1106 
1107  tAngent.resize(nb_gauss_pts, 3, false);
1108 
1109  switch (type) {
1110  case MBVERTEX:
1111  for (int dd = 0; dd != 3; dd++) {
1112  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1113  tAngent(gg, dd) = cblas_ddot(2, diff, 1, dofs[dd], 3);
1114  }
1115  }
1116  break;
1117  case MBEDGE:
1118  if (nb_dofs % 3) {
1119  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
1120  "Approximated field should be rank 3, i.e. vector in 3d space");
1121  }
1122  for (int dd = 0; dd != 3; dd++) {
1123  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1124  tAngent(gg, dd) +=
1125  cblas_ddot(nb_dofs / 3, &diff[gg * nb_approx_fun], 1, dofs[dd], 3);
1126  }
1127  }
1128  break;
1129  default:
1130  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
1131  "This operator can calculate tangent vector only on edge");
1132  }
1133 
1135 }
1136 
1138  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1140 
1141  if (type != MBEDGE)
1143 
1146  &tAngent[0], &tAngent[1], &tAngent[2]);
1147  const double l0 = t_m(i) * t_m(i);
1148 
1149  auto get_base_at_pts = [&](auto base) {
1151  &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
1152  &data.getN(base)(0, HVEC2));
1153  return t_h_curl;
1154  };
1155 
1156  auto get_tangent_at_pts = [&]() {
1158  &tangentAtGaussPt(0, 0), &tangentAtGaussPt(0, 1),
1159  &tangentAtGaussPt(0, 2));
1160  return t_m_at_pts;
1161  };
1162 
1163  auto calculate_squared_edge_length = [&]() {
1164  std::vector<double> l1;
1165  int nb_gauss_pts = tangentAtGaussPt.size1();
1166  if (nb_gauss_pts) {
1167  l1.resize(nb_gauss_pts);
1168  auto t_m_at_pts = get_tangent_at_pts();
1169  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1170  l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
1171  ++t_m_at_pts;
1172  }
1173  }
1174  return l1;
1175  };
1176 
1177  auto l1 = calculate_squared_edge_length();
1178 
1179  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
1180 
1181  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
1182  const size_t nb_gauss_pts = data.getN(base).size1();
1183  const size_t nb_dofs = data.getN(base).size2() / 3;
1184  if (nb_gauss_pts && nb_dofs) {
1185  auto t_h_curl = get_base_at_pts(base);
1186  int cc = 0;
1187  if (tangentAtGaussPt.size1() == nb_gauss_pts) {
1188  auto t_m_at_pts = get_tangent_at_pts();
1189  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1190  const double l0 = l1[gg];
1191  for (int ll = 0; ll != nb_dofs; ll++) {
1192  const double val = t_h_curl(0);
1193  const double a = val / l0;
1194  t_h_curl(i) = t_m_at_pts(i) * a;
1195  ++t_h_curl;
1196  ++cc;
1197  }
1198  ++t_m_at_pts;
1199  }
1200  } else {
1201  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1202  for (int ll = 0; ll != nb_dofs; ll++) {
1203  const double val = t_h_curl(0);
1204  const double a = val / l0;
1205  t_h_curl(i) = t_m(i) * a;
1206  ++t_h_curl;
1207  ++cc;
1208  }
1209  }
1210  }
1211 
1212  if (cc != nb_gauss_pts * nb_dofs)
1213  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Data inconsistency");
1214  }
1215  }
1216 
1218 }
1219 
1220 template <>
1221 template <>
1224  double *ptr = &*data.data().begin();
1225  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
1226 }
1227 
1228 template <>
1229 template <>
1231 OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(MatrixDouble &data) {
1232  double *ptr = &*data.data().begin();
1233  return FTensor::Tensor2<double *, 3, 3>(ptr, &ptr[1], &ptr[2], &ptr[3],
1234  &ptr[4], &ptr[5], &ptr[6], &ptr[7],
1235  &ptr[8], 9);
1236 }
1237 
1238 template <>
1240  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1242  if (data.getBase() == NOBASE)
1244  const unsigned int nb_gauss_pts = data.getN().size1();
1245  const unsigned int nb_base_functions = data.getN().size2();
1246  const unsigned int nb_dofs = data.getFieldData().size();
1247  if (!nb_dofs)
1249  auto t_n = data.getFTensor0N();
1250  auto t_val = getValAtGaussPtsTensor<3>(dataAtGaussPts);
1251  auto t_grad = getGradAtGaussPtsTensor<3, 3>(dataGradAtGaussPts);
1254  if (type == MBVERTEX &&
1255  data.getDiffN().data().size() == 3 * nb_base_functions) {
1256  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
1257  auto t_data = data.getFTensor1FieldData<3>();
1258  auto t_diff_n = data.getFTensor1DiffN<3>();
1259  unsigned int bb = 0;
1260  for (; bb != nb_dofs / 3; ++bb) {
1261  t_val(i) += t_data(i) * t_n;
1262  t_grad(i, j) += t_data(i) * t_diff_n(j);
1263  ++t_n;
1264  ++t_diff_n;
1265  ++t_data;
1266  }
1267  ++t_val;
1268  ++t_grad;
1269  for (; bb != nb_base_functions; ++bb) {
1270  ++t_n;
1271  }
1272  }
1273  } else {
1274  auto t_diff_n = data.getFTensor1DiffN<3>();
1275  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
1276  auto t_data = data.getFTensor1FieldData<3>();
1277  unsigned int bb = 0;
1278  for (; bb != nb_dofs / 3; ++bb) {
1279  t_val(i) += t_data(i) * t_n;
1280  t_grad(i, j) += t_data(i) * t_diff_n(j);
1281  ++t_n;
1282  ++t_diff_n;
1283  ++t_data;
1284  }
1285  ++t_val;
1286  ++t_grad;
1287  for (; bb != nb_base_functions; ++bb) {
1288  ++t_n;
1289  ++t_diff_n;
1290  }
1291  }
1292  }
1294 }
1295 
1296 template <>
1298  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1300  const unsigned int nb_gauss_pts = data.getN().size1();
1301  const unsigned int nb_base_functions = data.getN().size2();
1302  // bool constant_diff = false;
1303  const unsigned int nb_dofs = data.getFieldData().size();
1304  auto t_n = data.getFTensor0N();
1306  FTensor::Tensor0<double *>(&*dataAtGaussPts.data().begin(), 1);
1307  double *ptr = &*dataGradAtGaussPts.data().begin();
1308  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_grad(ptr, &ptr[1],
1309  &ptr[2]);
1311  if (type == MBVERTEX &&
1312  data.getDiffN().data().size() == 3 * nb_base_functions) {
1313  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
1314  auto t_data = data.getFTensor0FieldData();
1315  auto t_diff_n = data.getFTensor1DiffN<3>();
1316  unsigned int bb = 0;
1317  for (; bb != nb_dofs / 3; ++bb) {
1318  t_val += t_data * t_n;
1319  t_grad(i) += t_data * t_diff_n(i);
1320  ++t_n;
1321  ++t_diff_n;
1322  ++t_data;
1323  }
1324  ++t_val;
1325  ++t_grad;
1326  for (; bb != nb_base_functions; ++bb) {
1327  ++t_n;
1328  }
1329  }
1330  } else {
1331  auto t_diff_n = data.getFTensor1DiffN<3>();
1332  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
1333  auto t_data = data.getFTensor0FieldData();
1334  unsigned int bb = 0;
1335  for (; bb != nb_dofs / 3; ++bb) {
1336  t_val = t_data * t_n;
1337  t_grad(i) += t_data * t_diff_n(i);
1338  ++t_n;
1339  ++t_diff_n;
1340  ++t_data;
1341  }
1342  ++t_val;
1343  ++t_grad;
1344  for (; bb != nb_base_functions; ++bb) {
1345  ++t_n;
1346  ++t_diff_n;
1347  }
1348  }
1349  }
1351 }
1352 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
FTensor::Tensor2
Definition: Tensor2_value.hpp:17
MoFEM::OpSetContravariantPiolaTransform::j
FTensor::Index< 'j', 3 > j
Definition: DataOperators.hpp:218
MoFEM::OpSetHoContravariantPiolaTransform::i
FTensor::Index< 'i', 3 > i
Definition: DataOperators.hpp:242
MoFEM::OpSetHoInvJacH1::i
FTensor::Index< 'i', 3 > i
Definition: DataOperators.hpp:169
MoFEM::OpSetInvJacH1::i
FTensor::Index< 'i', 3 > i
Definition: DataOperators.hpp:124
MoFEM::DataForcesAndSourcesCore::EntData::getBase
FieldApproximationBase & getBase()
Get approximation base.
Definition: DataStructures.hpp:1281
MoFEM::OpGetCoordsAndNormalsOnPrism::tAngent2_at_GaussPtF4
MatrixDouble & tAngent2_at_GaussPtF4
Definition: DataOperators.hpp:480
MoFEM::OpSetHoCovariantPiolaTransform::hoInvJac
MatrixDouble & hoInvJac
Definition: DataOperators.hpp:259
LASTBASE
@ LASTBASE
Definition: definitions.h:161
MoFEM::OpSetHoContravariantPiolaTransform::k
FTensor::Index< 'k', 3 > k
Definition: DataOperators.hpp:244
HVEC0_2
@ HVEC0_2
Definition: definitions.h:267
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM::OpSetContravariantPiolaTransform::piolaDiffN
MatrixDouble piolaDiffN
Definition: DataOperators.hpp:228
ApproximationBaseNames
static const char *const ApproximationBaseNames[]
Definition: definitions.h:164
MoFEM::OpGetCoordsAndNormalsOnFace::tAngent2_at_GaussPt
MatrixDouble & tAngent2_at_GaussPt
Definition: DataOperators.hpp:450
MoFEM::OpSetCovariantPiolaTransformOnFace::tangent1AtGaussPt
const MatrixDouble & tangent1AtGaussPt
Definition: DataOperators.hpp:551
MoFEM::OpSetContravariantPiolaTransform::tJac
FTensor::Tensor2< double *, 3, 3 > tJac
Definition: DataOperators.hpp:216
MoFEM::OpSetInvJacHdivAndHcurl::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:270
MoFEM::DataOperator::getSymm
bool getSymm() const
Get if operator uses symmetry of DOFs or not.
Definition: DataOperators.hpp:97
NOBASE
@ NOBASE
Definition: definitions.h:151
MoFEM::Types::MatrixDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
MoFEM::OpSetCovariantPiolaTransform::i
FTensor::Index< 'i', 3 > i
Definition: DataOperators.hpp:288
MoFEM::OpSetInvJacHdivAndHcurl::k
FTensor::Index< 'k', 3 > k
Definition: DataOperators.hpp:148
MoFEM::OpSetHoInvJacHdivAndHcurl::j
FTensor::Index< 'j', 3 > j
Definition: DataOperators.hpp:189
cblas_dgemv
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
MoFEM::DataForcesAndSourcesCore::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: DataStructures.hpp:1220
HVEC0_1
@ HVEC0_1
Definition: definitions.h:264
CblasNoTrans
@ CblasNoTrans
Definition: cblas.h:11
MoFEM::OpSetCovariantPiolaTransform::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:380
FTensor::Tensor1
Definition: Tensor1_value.hpp:9
MoFEM::DataForcesAndSourcesCore::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: DataStructures.hpp:1300
MoFEM::OpSetInvJacHdivAndHcurl::diffHdivInvJac
MatrixDouble diffHdivInvJac
Definition: DataOperators.hpp:155
MoFEM::OpSetCovariantPiolaTransformOnFace::tAngent0
const VectorDouble & tAngent0
Definition: DataOperators.hpp:548
MoFEM::OpSetHoInvJacHdivAndHcurl::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:513
MoFEM::OpSetCovariantPiolaTransform::piolaN
MatrixDouble piolaN
Definition: DataOperators.hpp:297
MoFEM::OpSetHoCovariantPiolaTransform::piolaN
MatrixDouble piolaN
Definition: DataOperators.hpp:266
MoFEM::OpSetCovariantPiolaTransformOnFace::tAngent1
const VectorDouble & tAngent1
Definition: DataOperators.hpp:550
MoFEM::OpSetCovariantPiolaTransformOnFace::normalsAtGaussPts
const MatrixDouble & normalsAtGaussPts
Definition: DataOperators.hpp:547
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor0FieldData
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
Definition: DataStructures.cpp:402
HVEC0_0
@ HVEC0_0
Definition: definitions.h:261
MoFEM::OpSetCovariantPiolaTransformOnEdge::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:1137
MoFEM::DataForcesAndSourcesCore::EntData::getBBDiffNByOrderArray
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
Definition: DataStructures.cpp:849
HVEC1_1
@ HVEC1_1
Definition: definitions.h:265
MoFEM::determinantTensor3by3
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:905
HVEC1
@ HVEC1
Definition: definitions.h:255
FTensor::Number< 0 >
MoFEM::OpGetHoTangentOnEdge::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:1091
N0
static Number< 0 > N0
Definition: BasicFeTools.hpp:89
MoFEM::OpSetInvJacH1::diffNinvJac
MatrixDouble diffNinvJac
Definition: DataOperators.hpp:132
MoFEM::OpSetHoContravariantPiolaTransform::piolaDiffN
MatrixDouble piolaDiffN
Definition: DataOperators.hpp:250
MoFEM::OpSetHoInvJacHdivAndHcurl::i
FTensor::Index< 'i', 3 > i
Definition: DataOperators.hpp:188
MoFEM::OpSetCovariantPiolaTransform::tInvJac
FTensor::Tensor2< double *, 3, 3 > tInvJac
Definition: DataOperators.hpp:287
HVEC2_1
@ HVEC2_1
Definition: definitions.h:266
MoFEM::OpSetHoContravariantPiolaTransform::hoJac
MatrixDouble & hoJac
Definition: DataOperators.hpp:241
MoFEM::OpSetHoCovariantPiolaTransform::j
FTensor::Index< 'j', 3 > j
Definition: DataOperators.hpp:261
MoFEM::OpSetCovariantPiolaTransform::piolaDiffN
MatrixDouble piolaDiffN
Definition: DataOperators.hpp:298
m
static Index< 'm', 3 > m
Definition: BasicFeTools.hpp:77
MoFEM::OpSetHoInvJacHdivAndHcurl::k
FTensor::Index< 'k', 3 > k
Definition: DataOperators.hpp:190
MoFEM::OpSetContravariantPiolaTransformOnFace::normalShift
int normalShift
Shift in vector for linear geometry.
Definition: DataOperators.hpp:521
MoFEM::OpSetContravariantPiolaTransformOnFace::normalsAtGaussPtsRawPtr
const MatrixDouble * normalsAtGaussPtsRawPtr
Definition: DataOperators.hpp:510
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: DataStructures.cpp:587
MoFEM::OpSetInvJacHdivAndHcurl::i
FTensor::Index< 'i', 3 > i
Definition: DataOperators.hpp:146
HVEC2_2
@ HVEC2_2
Definition: definitions.h:269
HVEC1_2
@ HVEC1_2
Definition: definitions.h:268
MoFEM::OpSetHoContravariantPiolaTransform::piolaN
MatrixDouble piolaN
Definition: DataOperators.hpp:249
gm_rule.h
MoFEM::OpGetHoTangentOnEdge::tAngent
MatrixDouble & tAngent
Definition: DataOperators.hpp:572
MoFEM::OpSetHoInvJacH1::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:433
MoFEM::OpGetCoordsAndNormalsOnFace::cOords_at_GaussPt
MatrixDouble & cOords_at_GaussPt
Definition: DataOperators.hpp:447
MoFEM::OpSetInvJacHdivAndHcurl::tInvJac
FTensor::Tensor2< double *, 3, 3 > tInvJac
Definition: DataOperators.hpp:145
MoFEM::OpSetHoInvJacHdivAndHcurl::invHoJac
MatrixDouble & invHoJac
Definition: DataOperators.hpp:187
Spin
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:558
MoFEM::invertTensor3by3
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:886
MoFEM::OpGetCoordsAndNormalsOnPrism::nOrmals_at_GaussPtF3
MatrixDouble & nOrmals_at_GaussPtF3
Definition: DataOperators.hpp:474
MoFEM::OpGetDataAndGradient
Get field values and gradients at Gauss points.
Definition: DataOperators.hpp:307
convert.type
type
Definition: convert.py:66
MoFEM::OpSetInvJacHdivAndHcurl::j
FTensor::Index< 'j', 3 > j
Definition: DataOperators.hpp:147
MOFEM_IMPOSIBLE_CASE
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:127
CblasRowMajor
@ CblasRowMajor
Definition: cblas.h:10
FTensor::Index< 'i', 3 >
ElectroPhysiology::b
const double b
Definition: ElecPhysOperators.hpp:30
MoFEM::OpGetCoordsAndNormalsOnPrism::sPin
MatrixDouble sPin
Definition: DataOperators.hpp:496
MoFEM::OpSetCovariantPiolaTransformOnEdge::tangentAtGaussPt
const MatrixDouble & tangentAtGaussPt
Definition: DataOperators.hpp:586
lapack_wrap.h
MoFEM::DataOperator::DataOperator
DataOperator(const bool symm=true)
Definition: DataOperators.cpp:33
MoFEM::OpGetCoordsAndNormalsOnFace::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:674
MoFEM::OpGetCoordsAndNormalsOnFace::nOrmals_at_GaussPt
MatrixDouble & nOrmals_at_GaussPt
Definition: DataOperators.hpp:448
MoFEM::OpSetContravariantPiolaTransform::piolaN
MatrixDouble piolaN
Definition: DataOperators.hpp:227
MoFEM::OpSetCovariantPiolaTransform::k
FTensor::Index< 'k', 3 > k
Definition: DataOperators.hpp:290
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: DataStructures.hpp:1459
MoFEM::OpSetInvJacH1::j
FTensor::Index< 'j', 3 > j
Definition: DataOperators.hpp:125
N1
static Number< 1 > N1
Definition: BasicFeTools.hpp:90
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Definition: DataStructures.cpp:607
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEM::OpSetInvJacH1::tInvJac
FTensor::Tensor2< double *, 3, 3 > tInvJac
Definition: DataOperators.hpp:123
FTensor::Tensor0
Definition: Tensor0.hpp:17
MoFEM::OpSetHoInvJacH1::j
FTensor::Index< 'j', 3 > j
Definition: DataOperators.hpp:170
MoFEM::OpSetHoContravariantPiolaTransform::detHoJac
VectorDouble & detHoJac
Definition: DataOperators.hpp:240
cblas.h
MoFEM::OpSetCovariantPiolaTransform::j
FTensor::Index< 'j', 3 > j
Definition: DataOperators.hpp:289
MoFEM::OpGetCoordsAndNormalsOnPrism::tAngent1_at_GaussPtF4
MatrixDouble & tAngent1_at_GaussPtF4
Definition: DataOperators.hpp:479
MoFEM::DataOperator::doWork
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: DataOperators.hpp:45
MoFEM::DataOperator::opLhs
virtual MoFEMErrorCode opLhs(DataForcesAndSourcesCore &row_data, DataForcesAndSourcesCore &col_data)
Definition: DataOperators.cpp:111
MoFEM::OpSetContravariantPiolaTransform::i
FTensor::Index< 'i', 3 > i
Definition: DataOperators.hpp:217
MoFEM::DataForcesAndSourcesCore
data structure for finite element entity
Definition: DataStructures.hpp:52
MoFEM::OpSetCovariantPiolaTransformOnFace::tangent0AtGaussPt
const MatrixDouble & tangent0AtGaussPt
Definition: DataOperators.hpp:549
MoFEM::OpSetHoCovariantPiolaTransform::k
FTensor::Index< 'k', 3 > k
Definition: DataOperators.hpp:262
MoFEM::DataOperator::opRhs
virtual MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data, const bool error_if_no_base=false)
Definition: DataOperators.cpp:171
MoFEM::OpGetCoordsAndNormalsOnPrism::nOrmals_at_GaussPtF4
MatrixDouble & nOrmals_at_GaussPtF4
Definition: DataOperators.hpp:478
MoFEM::OpSetContravariantPiolaTransformOnFace::normalRawPtr
const VectorDouble * normalRawPtr
Definition: DataOperators.hpp:509
MoFEM::getFTensor0FromVec
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:143
MoFEM::OpGetCoordsAndNormalsOnPrism::tAngent2_at_GaussPtF3
MatrixDouble & tAngent2_at_GaussPtF3
Definition: DataOperators.hpp:476
MoFEM::OpGetCoordsAndNormalsOnPrism::tAngent1_at_GaussPtF3
MatrixDouble & tAngent1_at_GaussPtF3
Definition: DataOperators.hpp:475
MoFEM::OpGetCoordsAndNormalsOnFace::calculateNormals
MoFEMErrorCode calculateNormals()
Definition: DataOperators.cpp:759
MoFEM::OpSetCovariantPiolaTransformOnFace::nOrmal
const VectorDouble & nOrmal
Definition: DataOperators.hpp:546
MoFEM::DataForcesAndSourcesCore::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: DataStructures.hpp:1040
MoFEM::OpSetHoContravariantPiolaTransform::j
FTensor::Index< 'j', 3 > j
Definition: DataOperators.hpp:243
MoFEM::DataForcesAndSourcesCore::EntData::getBBDiffNMap
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
Definition: DataStructures.cpp:804
MoFEM::OpSetContravariantPiolaTransform::vOlume
double & vOlume
Definition: DataOperators.hpp:214
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::OpSetContravariantPiolaTransform::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:314
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
MoFEM::OpSetHoInvJacH1::invHoJac
MatrixDouble & invHoJac
Definition: DataOperators.hpp:168
MoFEM::OpGetCoordsAndNormalsOnPrism::cOords_at_GaussPtF4
MatrixDouble & cOords_at_GaussPtF4
Definition: DataOperators.hpp:477
MoFEM::OpSetCovariantPiolaTransformOnFace::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:988
MoFEM::OpSetCovariantPiolaTransformOnEdge::tAngent
const VectorDouble & tAngent
Definition: DataOperators.hpp:585
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
ContactOps::normal
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:75
MoFEM::OpSetContravariantPiolaTransform::k
FTensor::Index< 'k', 3 > k
Definition: DataOperators.hpp:219
MoFEM::DataForcesAndSourcesCore::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: DataStructures.hpp:1288
MoFEM::OpGetCoordsAndNormalsOnPrism::calculateNormals
MoFEMErrorCode calculateNormals()
Definition: DataOperators.cpp:871
MoFEM::OpSetHoCovariantPiolaTransform::piolaDiffN
MatrixDouble piolaDiffN
Definition: DataOperators.hpp:267
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:150
MoFEM::OpSetHoContravariantPiolaTransform::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:561
HVEC2_0
@ HVEC2_0
Definition: definitions.h:263
MoFEM::OpSetHoCovariantPiolaTransform::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:618
MoFEM::DataForcesAndSourcesCore::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: DataStructures.hpp:1235
MoFEM::OpGetCoordsAndNormalsOnPrism::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:787
cblas_ddot
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
ElectroPhysiology::c
const double c
Definition: ElecPhysOperators.hpp:31
MoFEM::OpGetCoordsAndNormalsOnFace::tAngent1_at_GaussPt
MatrixDouble & tAngent1_at_GaussPt
Definition: DataOperators.hpp:449
MoFEM::OpSetInvJacH1::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:200
MoFEM::DataForcesAndSourcesCore::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: DataStructures.hpp:60
MoFEM::OpSetHoInvJacH1::diffNinvJac
MatrixDouble diffNinvJac
Definition: DataOperators.hpp:173
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
HVEC2
@ HVEC2
Definition: definitions.h:255
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DataOperator::doEntities
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Definition: DataOperators.hpp:75
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor1FieldData
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
Definition: DataStructures.hpp:1256
MoFEM::Types::VectorDouble
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MoFEM::OpSetHoCovariantPiolaTransform::i
FTensor::Index< 'i', 3 > i
Definition: DataOperators.hpp:260
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
MoFEM::OpSetContravariantPiolaTransformOnFace::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: DataOperators.cpp:896
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:128
MoFEM::OpSetHoInvJacHdivAndHcurl::diffHdivInvJac
MatrixDouble diffHdivInvJac
Definition: DataOperators.hpp:194
I
static Index< 'I', 3 > I
Definition: BasicFeTools.hpp:82
HVEC1_0
@ HVEC1_0
Definition: definitions.h:262
FTensor::dd
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
HVEC0
@ HVEC0
Definition: definitions.h:255
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor1DiffN
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: DataStructures.cpp:415
MoFEM::OpGetCoordsAndNormalsOnPrism::cOords_at_GaussPtF3
MatrixDouble & cOords_at_GaussPtF3
Definition: DataOperators.hpp:473