v0.14.0
EntitiesFieldData.cpp
Go to the documentation of this file.
1 /** \file EntitiesFieldData.cpp
2 \brief Implementation for Data Structures in Forces and Sources
3 
4 */
5 
6 
7 
8 namespace MoFEM {
9 
10 EntitiesFieldData::EntData::EntData(const bool allocate_base_matrices)
11  : sEnse(0), oRder(0), bAse(NOBASE), entDataBitRefLevel(),
12  N(baseFunctionsAndBaseDerivatives[ZeroDerivative]),
13  diffN(baseFunctionsAndBaseDerivatives[FirstDerivative]) {
14  if (allocate_base_matrices) {
15 
16  for (auto d = 0; d != LastDerivative; ++d) {
17  for (int b = 0; b != LASTBASE; ++b) {
19  }
20  }
21  }
22 }
23 
24 int EntitiesFieldData::EntData::getSense() const { return sEnse; }
25 
26 boost::shared_ptr<MatrixDouble> &
28  const BaseDerivatives direvatie) {
29  return baseFunctionsAndBaseDerivatives[direvatie][base];
30 }
31 
32 boost::shared_ptr<MatrixDouble> &
34  return N[base];
35 }
36 
37 boost::shared_ptr<MatrixDouble> &EntitiesFieldData::EntData::getDiffNSharedPtr(
38  const FieldApproximationBase base) {
39  return diffN[base];
40 }
41 
42 static void constructor_data(EntitiesFieldData *data, const EntityType type) {
43 
45 
46  auto set_default = [&]() {
47  std::array<size_t, MBMAXTYPE> count;
48  std::fill(count.begin(), count.end(), 0);
49  const int dim_type = moab::CN::Dimension(type);
50  data->dataOnEntities[MBVERTEX].resize(1);
51  if (type != MBVERTEX) {
52  for (auto dd = dim_type; dd > 0; --dd) {
53  int nb_ents = moab::CN::NumSubEntities(type, dd);
54  for (int ii = 0; ii != nb_ents; ++ii) {
55  auto sub_ent_type = moab::CN::SubEntityType(type, dd, ii);
56  count[sub_ent_type] = nb_ents;
57  }
58  for (auto tt = moab::CN::TypeDimensionMap[dd].first;
59  tt <= moab::CN::TypeDimensionMap[dd].second; ++tt) {
60  data->dataOnEntities[tt].resize(count[tt]);
61  }
62  }
63  }
64  };
65 
66  switch (type) {
67  case MBENTITYSET:
68  break;
69 
70  default:
71  set_default();
72  }
73 }
74 
76  constructor_data(this, type);
77 }
78 
81  constructor_data(this, type);
83 }
84 
85 static void
87  const boost::shared_ptr<EntitiesFieldData> &data_ptr) {
88 
90  using DerivedEntData = DerivedEntitiesFieldData::DerivedEntData;
91 
92  for (int tt = MBVERTEX; tt != MBMAXTYPE; ++tt) {
93  auto &ent_data = data_ptr->dataOnEntities[tt];
94  auto &derived_ent_data = derived_data->dataOnEntities[tt];
95  for (auto c = derived_ent_data.size(); c < ent_data.size(); ++c) {
96  boost::shared_ptr<EntData> ent_data_ptr(data_ptr, &ent_data[c]);
97  derived_ent_data.push_back(new DerivedEntData(ent_data_ptr));
98  }
99  derived_ent_data.resize(ent_data.size());
100  }
101 }
102 
104  const boost::shared_ptr<EntitiesFieldData> &data_ptr)
105  : EntitiesFieldData(), dataPtr(data_ptr) {
107 }
108 
113 }
114 
117  sPace = NOSPACE;
118  bAse = NOBASE;
119  fieldEntities.resize(0, false);
120  iNdices.resize(0, false);
121  localIndices.resize(0, false);
122  dOfs.resize(0, false);
123  fieldData.resize(0, false);
125 }
126 
129  for (EntityType t = MBVERTEX; t != MBMAXTYPE; t++)
130  for (auto &e : dataOnEntities[t])
131  CHKERR e.resetFieldDependentData();
133 }
134 
137  const FieldApproximationBase base) {
139  auto make_swap = [](boost::shared_ptr<MatrixDouble> &ptr,
140  boost::shared_ptr<MatrixDouble> &ptrBB,
141  boost::shared_ptr<MatrixDouble> &swap_ptr) {
142  if (swap_ptr) {
143  ptr = swap_ptr;
144  swap_ptr.reset();
145  } else {
146  swap_ptr = ptr;
147  ptr = ptrBB;
148  }
149  };
150  make_swap(getNSharedPtr(base), getBBNSharedPtr(field_name), swapBaseNPtr);
151  make_swap(getDiffNSharedPtr(base), getBBDiffNSharedPtr(field_name),
152  swapBaseDiffNPtr);
154 }
155 
157  const FieldApproximationBase base) {
159  // Note: Do not swap bases on entities sets
160  for (int tt = MBVERTEX; tt != MBENTITYSET; ++tt) {
161  auto &ent_data = dataOnEntities[tt];
162  for (auto &side_data : ent_data)
163  CHKERR side_data.baseSwap(field_name, base);
164  }
166 }
167 
169  const std::string &field_name, const FieldApproximationBase base) {
171  auto make_swap = [](boost::shared_ptr<MatrixDouble> &ptr,
172  boost::shared_ptr<MatrixDouble> &ptrBB,
173  boost::shared_ptr<MatrixDouble> &swap_ptr) {
174  if (swap_ptr) {
175  ptr = swap_ptr;
176  swap_ptr.reset();
177  } else {
178  swap_ptr = ptr;
179  ptr = ptrBB;
180  }
181  };
182  make_swap(getDerivedNSharedPtr(base), getBBNSharedPtr(field_name),
183  swapBaseNPtr);
184  make_swap(getDerivedDiffNSharedPtr(base), getBBDiffNSharedPtr(field_name),
185  swapBaseDiffNPtr);
187 }
188 
190  const boost::shared_ptr<EntitiesFieldData::EntData> &ent_data_ptr)
191  : EntitiesFieldData::EntData(false), entDataPtr(ent_data_ptr) {}
192 
194  return entDataPtr->getSense();
195 }
196 
197 boost::shared_ptr<MatrixDouble> &
199  const FieldApproximationBase base, const BaseDerivatives derivative) {
200  if (baseFunctionsAndBaseDerivatives[derivative][base])
201  return baseFunctionsAndBaseDerivatives[derivative][base];
202  else
203  return entDataPtr->getNSharedPtr(base, derivative);
204 }
205 
206 boost::shared_ptr<MatrixDouble> &
208  const FieldApproximationBase base) {
209  if (N[base])
210  return N[base];
211  else
212  return entDataPtr->getNSharedPtr(base);
213 }
214 
215 boost::shared_ptr<MatrixDouble> &
217  const FieldApproximationBase base) {
218  if (diffN[base])
219  return diffN[base];
220  else
221  return entDataPtr->getDiffNSharedPtr(base);
222 }
223 const boost::shared_ptr<MatrixDouble> &
225  const FieldApproximationBase base) const {
226  if (N[base])
227  return N[base];
228  else
229  return entDataPtr->getNSharedPtr(base);
230 }
231 const boost::shared_ptr<MatrixDouble> &
233  const FieldApproximationBase base) const {
234  if (diffN[base])
235  return diffN[base];
236  else
237  return entDataPtr->getDiffNSharedPtr(base);
238 }
239 
240 std::ostream &operator<<(std::ostream &os,
241  const EntitiesFieldData::EntData &e) {
242  os << "sEnse: " << e.getSense() << std::endl
243  << "oRder: " << e.getOrder() << std::endl
244  << "global indices: " << e.getIndices() << std::endl
245  << "local indices: " << e.getLocalIndices() << std::endl;
246  // FIXME: precision should not be set here
247  os << "fieldData: " << std::fixed << std::setprecision(2) << e.getFieldData()
248  << std::endl;
249  MatrixDouble base = const_cast<EntitiesFieldData::EntData &>(e).getN();
250  MatrixDouble diff_base =
251  const_cast<EntitiesFieldData::EntData &>(e).getDiffN();
252  const double eps = 1e-6;
253  for (unsigned int ii = 0; ii != base.size1(); ii++) {
254  for (unsigned int jj = 0; jj != base.size2(); jj++) {
255  if (fabs(base(ii, jj)) < eps)
256  base(ii, jj) = 0;
257  }
258  }
259  for (unsigned int ii = 0; ii != diff_base.size1(); ii++) {
260  for (unsigned int jj = 0; jj != diff_base.size2(); jj++) {
261  if (fabs(diff_base(ii, jj)) < eps)
262  diff_base(ii, jj) = 0;
263  }
264  }
265  os << "N: " << std::fixed << base << std::endl
266  << "diffN: " << std::fixed << diff_base;
267  return os;
268 }
269 
270 std::ostream &operator<<(std::ostream &os, const EntitiesFieldData &e) {
271  for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
272  for (unsigned int nn = 0; nn < e.dataOnEntities[t].size(); nn++) {
273  os << "dataOnEntities[" << moab::CN::EntityTypeName(t) << "][" << nn
274  << "]" << std::endl
275  << e.dataOnEntities[t][nn] << std::endl;
276  }
277  }
278  return os;
279 }
280 
281 /** \name Specializations for H1/L2 */
282 
283 /**@{*/
284 
285 template <>
287 EntitiesFieldData::EntData::getFTensor1FieldData<3>() {
288 #ifndef NDEBUG
289  for (auto &d : dOfs) {
290  if (d) {
291  if (d->getNbOfCoeffs() != 3) {
292  std::stringstream s;
293  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
294  s << " but you ask for tensor rank 1 dimension 3";
296  }
297  break;
298  }
299  }
300 #endif
301  double *ptr = &*fieldData.data().begin();
302  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
303  &ptr[2]);
304 }
305 
306 template <>
308 EntitiesFieldData::EntData::getFTensor1FieldData<2>() {
309 #ifndef NDEBUG
310  for (auto &d : dOfs) {
311  if (d) {
312  if (d->getNbOfCoeffs() != 2) {
313  std::stringstream s;
314  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
315  s << " but you ask for tensor rank 1 dimension 2";
317  }
318  break;
319  }
320  }
321 #endif
322  double *ptr = &*fieldData.data().begin();
323  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(ptr, &ptr[1]);
324 }
325 
326 template <>
328 EntitiesFieldData::EntData::getFTensor1FieldData<1>() {
329 #ifndef NDEBUG
330  for (auto &d : dOfs) {
331  if (d) {
332  if (d->getNbOfCoeffs() != 1) {
333  std::stringstream s;
334  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
335  s << " but you ask for tensor rank 1 dimension 1";
337  }
338  break;
339  }
340  }
341 #endif
342  double *ptr = &*fieldData.data().begin();
344 }
345 
346 template <>
348 EntitiesFieldData::EntData::getFTensor2FieldData<1, 1>() {
349 #ifndef NDEBUG
350  for (auto &d : dOfs) {
351  if (d) {
352  if (d->getNbOfCoeffs() != 1) {
353  std::stringstream s;
354  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
355  s << " but you ask for tensor rank 2 dimensions 1 by 1 so 1 "
356  "coefficients "
357  "is expected";
359  }
360  break;
361  }
362  }
363 #endif
364  double *ptr = &*fieldData.data().begin();
366 }
367 
368 template <>
370 EntitiesFieldData::EntData::getFTensor2FieldData<1, 2>() {
371 #ifndef NDEBUG
372  for (auto &d : dOfs) {
373  if (d) {
374  if (d->getNbOfCoeffs() != 2) {
375  std::stringstream s;
376  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
377  s << " but you ask for tensor rank 2 dimensions 1 by 2 so 2 "
378  "coefficients "
379  "is expected";
381  }
382  break;
383  }
384  }
385 #endif
386  double *ptr = &*fieldData.data().begin();
387  return FTensor::Tensor2<FTensor::PackPtr<double *, 2>, 1, 2>(ptr, &ptr[1]);
388 }
389 
390 template <>
392 EntitiesFieldData::EntData::getFTensor2FieldData<1, 3>() {
393 #ifndef NDEBUG
394  for (auto &d : dOfs) {
395  if (d) {
396  if (d->getNbOfCoeffs() != 3) {
397  std::stringstream s;
398  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
399  s << " but you ask for tensor rank 2 dimensions 1 by 3 so 3 "
400  "coefficients "
401  "is expected";
403  }
404  break;
405  }
406  }
407 #endif
408  double *ptr = &*fieldData.data().begin();
409  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 1, 3>(ptr, &ptr[1],
410  &ptr[2]);
411 }
412 
413 template <>
415 EntitiesFieldData::EntData::getFTensor2FieldData<2, 2>() {
416 #ifndef NDEBUG
417  for(auto &d : dOfs) {
418  if(d) {
419  if(d->getNbOfCoeffs() != 4) {
420  std::stringstream s;
421  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
422  s << " but you ask for tensor rank 2 dimensions 2 by 2 so 4 coefficients "
423  "is expected";
425  }
426  break;
427  }
428  }
429 #endif
430  double *ptr = &*fieldData.data().begin();
432  ptr, &ptr[1], &ptr[2], &ptr[3]);
433 }
434 
435 template <>
437 EntitiesFieldData::EntData::getFTensor2FieldData<3, 3>() {
438 #ifndef NDEBUG
439  for (auto &d : dOfs) {
440  if (d) {
441  if (d->getNbOfCoeffs() != 9) {
442  std::stringstream s;
443  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
444  s << " but you ask for tensor rank 2 dimensions 3 by 3 so 9 "
445  "coefficients "
446  "is expected";
448  }
449  break;
450  }
451  }
452 #endif
453  double *ptr = &*fieldData.data().begin();
455  ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
456  &ptr[8]);
457 }
458 
459 template <>
461 EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<3>() {
462 #ifndef NDEBUG
463  for (auto &d : dOfs) {
464  if (d) {
465  if (d->getNbOfCoeffs() != 6) {
466  std::stringstream s;
467  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
468  s << " but you ask for symmetric tensor rank 2 dimensions 3 by 3 so 6 "
469  "coefficients "
470  "is expected";
472  }
473  break;
474  }
475  }
476 #endif
477  double *ptr = &*fieldData.data().begin();
479  ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
480 }
481 
482 template <>
484 EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<2>() {
485 #ifndef NDEBUG
486  for (auto &d : dOfs) {
487  if (d) {
488  if (d->getNbOfCoeffs() != 3) {
489  std::stringstream s;
490  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
491  s << " but you ask for symmetric tensor rank 2 dimensions 2 by 2 so 3 "
492  "coefficients "
493  "is expected";
495  }
496  break;
497  }
498  }
499 #endif
500  double *ptr = &*fieldData.data().begin();
502  ptr, &ptr[1], &ptr[2]);
503 }
504 
507 #ifndef NDEBUG
508  for (auto &d : dOfs) {
509  if (d) {
510  if (d->getNbOfCoeffs() != 1) {
511  std::stringstream s;
512  s << "Wrong number of coefficients is " << d->getNbOfCoeffs();
513  s << " but expected scalar field, tensor of rank 0";
515  }
516  break;
517  }
518  }
519 #endif
521  &*fieldData.data().begin());
522 }
523 
524 template <int Tensor_Dim>
527  const FieldApproximationBase base) {
528  std::stringstream s;
529  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
530  THROW_MESSAGE(s.str());
532 }
533 
534 template <int Tensor_Dim>
537  return getFTensor1DiffN<Tensor_Dim>(bAse);
538 }
539 
540 /**
541  * \brief Get spatial derivative of base function tensor for dimension 3d
542  */
543 template <>
545 EntitiesFieldData::EntData::getFTensor1DiffN<3>(
546  const FieldApproximationBase base) {
547  double *ptr = &*getDiffN(base).data().begin();
548  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
549  &ptr[2]);
550 }
551 
552 /**
553  * \brief Get spatial derivative of base function tensor for dimension 3d
554  */
555 template <>
557 EntitiesFieldData::EntData::getFTensor1DiffN<3>() {
558  return getFTensor1DiffN<3>(bAse);
559 }
560 
561 /**
562  * \brief Get spatial derivative of base function tensor for dimension 2d
563  */
564 template <>
566 EntitiesFieldData::EntData::getFTensor1DiffN<2>(
567  const FieldApproximationBase base) {
568  double *ptr = &*getDiffN(base).data().begin();
569  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(ptr, &ptr[1]);
570 }
571 
572 /**
573  * \brief Get spatial derivative of base function tensor for dimension 2d
574  */
575 template <>
577 EntitiesFieldData::EntData::getFTensor1DiffN<2>() {
578  return getFTensor1DiffN<2>(bAse);
579 }
580 
581 template <int Tensor_Dim>
584  const int gg, const int bb) {
585  std::stringstream s;
586  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
587  THROW_MESSAGE(s.str());
589 }
590 
591 /**
592  * \brief Get spatial derivative of base function tensor for dimension 3d
593  */
594 template <>
596 EntitiesFieldData::EntData::getFTensor1DiffN<3>(
597  const FieldApproximationBase base, const int gg, const int bb) {
598  double *ptr = &getDiffN(base)(gg, 3 * bb);
599  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
600  &ptr[2]);
601 }
602 
603 /**
604  * \brief Get spatial derivative of base function tensor for dimension 3d
605  */
606 template <>
608 EntitiesFieldData::EntData::getFTensor1DiffN<3>(const int gg, const int bb) {
609  return getFTensor1DiffN<3>(bAse, gg, bb);
610 }
611 
612 /**
613  * \brief Get spatial derivative of base function tensor for dimension 2d
614  */
615 template <>
617 EntitiesFieldData::EntData::getFTensor1DiffN<2>(
618  const FieldApproximationBase base, const int gg, const int bb) {
619  double *ptr = &getDiffN(base)(gg, 2 * bb);
620  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(ptr, &ptr[1]);
621 }
622 
623 /**
624  * \brief Get spatial derivative of base function tensor for dimension 2d
625  */
626 template <>
628 EntitiesFieldData::EntData::getFTensor1DiffN<2>(const int gg, const int bb) {
629  return getFTensor1DiffN<2>(bAse, gg, bb);
630 }
631 
632 /**@}*/
633 
634 /** \name Specializations for HDiv/HCrul */
635 
636 /**@{*/
637 
638 template <int Tensor_Dim>
641  std::stringstream s;
642  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
643  THROW_MESSAGE(s.str());
645 }
646 
647 template <int Tensor_Dim>
650  const int gg, const int bb) {
651  std::stringstream s;
652  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
653  THROW_MESSAGE(s.str());
655 }
656 
657 template <int Tensor_Dim0, int Tensor_Dim1>
659  Tensor_Dim0, Tensor_Dim1>
661  std::stringstream s;
662  s << "Template for tensor dimension " << Tensor_Dim0 << "x" << Tensor_Dim1
663  << " not implemented";
664  THROW_MESSAGE(s.str());
666 }
667 
668 template <int Tensor_Dim0, int Tensor_Dim1>
670  Tensor_Dim0, Tensor_Dim1>
672  const int gg, const int bb) {
673  std::stringstream s;
674  s << "Template for tensor dimension " << Tensor_Dim0 << "x" << Tensor_Dim1
675  << " not implemented";
676  THROW_MESSAGE(s.str());
678 }
679 
680 template <>
682 EntitiesFieldData::EntData::getFTensor1N<3>(FieldApproximationBase base) {
683  double *t_n_ptr = &*getN(base).data().begin();
684  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(t_n_ptr, // HVEC0
685  &t_n_ptr[HVEC1],
686  &t_n_ptr[HVEC2]);
687 }
688 
689 template <>
691 EntitiesFieldData::EntData::getFTensor1N<3>(FieldApproximationBase base,
692  const int gg, const int bb) {
693  double *t_n_ptr = &getN(base)(gg, 3 * bb);
694  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(t_n_ptr, // HVEC0
695  &t_n_ptr[HVEC1],
696  &t_n_ptr[HVEC2]);
697 }
698 
699 template <>
701 EntitiesFieldData::EntData::getFTensor2DiffN<3, 3>(
702  FieldApproximationBase base) {
703  double *t_diff_n_ptr = &*getDiffN(base).data().begin();
705  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC0_2],
706  &t_diff_n_ptr[HVEC1_0], &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC1_2],
707  &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1], &t_diff_n_ptr[HVEC2_2]);
708 }
709 
710 template <>
712 EntitiesFieldData::EntData::getFTensor2DiffN<3, 3>(FieldApproximationBase base,
713  const int gg, const int bb) {
714  double *t_diff_n_ptr = &getDiffN(base)(gg, 9 * bb);
716  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC0_2],
717  &t_diff_n_ptr[HVEC1_0], &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC1_2],
718  &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1], &t_diff_n_ptr[HVEC2_2]);
719 }
720 
721 template <>
723 EntitiesFieldData::EntData::getFTensor2DiffN<3, 2>(
724  FieldApproximationBase base) {
725  double *t_diff_n_ptr = &*getDiffN(base).data().begin();
727  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC1_0],
728  &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1]);
729 }
730 
731 template <>
733 EntitiesFieldData::EntData::getFTensor2DiffN<3, 2>(FieldApproximationBase base,
734  const int gg, const int bb) {
735  double *t_diff_n_ptr = &getDiffN(base)(gg, 6 * bb);
737  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC1_0],
738  &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1]);
739 }
740 
741 template <>
744  double *ptr = &(getN(base, BaseDerivatives::SecondDerivative))(0, 0);
746 
747  &ptr[2 * HVEC0_0 + 0], &ptr[2 * HVEC0_0 + 1], &ptr[2 * HVEC0_1 + 0],
748  &ptr[2 * HVEC0_1 + 1],
749 
750  &ptr[2 * HVEC1_0 + 0], &ptr[2 * HVEC1_0 + 1], &ptr[2 * HVEC1_1 + 0],
751  &ptr[2 * HVEC1_1 + 1],
752 
753  &ptr[2 * HVEC2_0 + 0], &ptr[2 * HVEC2_0 + 1], &ptr[2 * HVEC2_1 + 0],
754  &ptr[2 * HVEC2_1 + 1]
755 
756  };
757 }
758 
759 template <int Tensor_Dim0, int Tensor_Dim1>
761  Tensor_Dim0, Tensor_Dim1>
763  std::stringstream s;
764  s << "Template for tensor dimension " << Tensor_Dim0 << ", " << Tensor_Dim1
765  << " not implemented";
766  THROW_MESSAGE(s.str());
768  Tensor_Dim0, Tensor_Dim1>();
769 }
770 
771 template <>
773 EntitiesFieldData::EntData::getFTensor2N<3, 3>(FieldApproximationBase base) {
774  double *t_n_ptr = &*(getN(base).data().begin());
776 
777  &t_n_ptr[HVEC0], &t_n_ptr[HVEC1], &t_n_ptr[HVEC2],
778 
779  &t_n_ptr[3 + HVEC0], &t_n_ptr[3 + HVEC1], &t_n_ptr[3 + HVEC2],
780 
781  &t_n_ptr[6 + HVEC0], &t_n_ptr[6 + HVEC1], &t_n_ptr[6 + HVEC2]
782 
783  );
784 }
785 
786 template <int Tensor_Dim0, int Tensor_Dim1>
788  Tensor_Dim0, Tensor_Dim1>
790  const int gg, const int bb) {
791  std::stringstream s;
792  s << "Template for tensor dimension " << Tensor_Dim0 << ", " << Tensor_Dim1
793  << " not implemented";
794  THROW_MESSAGE(s.str());
796  Tensor_Dim0, Tensor_Dim1>();
797 }
798 
799 template <>
801 EntitiesFieldData::EntData::getFTensor2N<3, 3>(FieldApproximationBase base,
802  const int gg, const int bb) {
803  double *t_n_ptr = &getN(base)(gg, 9 * bb);
805 
806  &t_n_ptr[HVEC0], &t_n_ptr[HVEC1], &t_n_ptr[HVEC2],
807 
808  &t_n_ptr[3 + HVEC0], &t_n_ptr[3 + HVEC1], &t_n_ptr[3 + HVEC2],
809 
810  &t_n_ptr[6 + HVEC0], &t_n_ptr[6 + HVEC1], &t_n_ptr[6 + HVEC2]
811 
812  );
813 }
814 
815 /**@}*/
816 
817 /** \name Bernstein-Bezier base only functions */
818 
819 /**@{*/
820 
821 boost::shared_ptr<MatrixInt> &
823  const std::string &field_name) {
824  return bbAlphaIndices[field_name];
825 }
826 
827 boost::shared_ptr<MatrixDouble> &
829  return bbN[field_name];
830 }
831 
832 /**
833  * Get shared pointer to BB base base functions
834  */
835 const boost::shared_ptr<MatrixDouble> &
837  const std::string &field_name) const {
838  return bbN.at(field_name);
839 }
840 
841 /**
842  * Get shared pointer to BB derivatives of base base functions
843  */
844 boost::shared_ptr<MatrixDouble> &
846  return bbDiffN[field_name];
847 }
848 
849 /**
850  * Get shared pointer to derivatives of BB base base functions
851  */
852 const boost::shared_ptr<MatrixDouble> &
854  const std::string &field_name) const {
855  return bbDiffN.at(field_name);
856 }
857 
858 std::map<std::string, boost::shared_ptr<MatrixInt>> &
860  return bbAlphaIndices;
861 }
862 
863 std::map<std::string, boost::shared_ptr<MatrixDouble>> &
865  return bbN;
866 }
867 
868 std::map<std::string, boost::shared_ptr<MatrixDouble>> &
870  return bbDiffN;
871 }
872 
873 boost::shared_ptr<MatrixInt> &
875  return bbAlphaIndicesByOrder[o];
876 }
877 
878 boost::shared_ptr<MatrixDouble> &
880  return bbNByOrder[o];
881 }
882 
883 boost::shared_ptr<MatrixDouble> &
885  return bbDiffNByOrder[o];
886 }
887 
888 std::array<boost::shared_ptr<MatrixInt>,
891  return bbAlphaIndicesByOrder;
892 }
893 
894 std::array<boost::shared_ptr<MatrixDouble>,
897  return bbNByOrder;
898 }
899 
900 std::array<boost::shared_ptr<MatrixDouble>,
903  return bbDiffNByOrder;
904 }
905 
906 boost::shared_ptr<MatrixInt> &
908  const std::string &field_name) {
909  return entDataPtr->getBBAlphaIndicesSharedPtr(field_name);
910 }
911 
912 boost::shared_ptr<MatrixDouble> &
914  const std::string &field_name) {
915  return entDataPtr->getBBNSharedPtr(field_name);
916 }
917 
918 const boost::shared_ptr<MatrixDouble> &
920  const std::string &field_name) const {
921  return entDataPtr->getBBNSharedPtr(field_name);
922 }
923 
924 /**
925  * Get shared pointer to BB derivatives of base base functions
926  */
927 boost::shared_ptr<MatrixDouble> &
929  const std::string &field_name) {
930  return entDataPtr->getBBDiffNSharedPtr(field_name);
931 }
932 
933 /**
934  * Get shared pointer to derivatives of BB base base functions
935  */
936 const boost::shared_ptr<MatrixDouble> &
938  const std::string &field_name) const {
939  return entDataPtr->getBBDiffNSharedPtr(field_name);
940 }
941 
942 /**@}*/
943 
945  return entDataBitRefLevel;
946 }
947 
948 std::vector<BitRefLevel> &
950  return entDataPtr->getEntDataBitRefLevel();
951 }
952 
953 } // namespace MoFEM
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::EntitiesFieldData::EntData::resetFieldDependentData
MoFEMErrorCode resetFieldDependentData()
Definition: EntitiesFieldData.cpp:115
MoFEM::EntitiesFieldData::EntData::getBBAlphaIndicesSharedPtr
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
Definition: EntitiesFieldData.cpp:822
MoFEM::EntitiesFieldData::EntData::getBBDiffNMap
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
get hash map of direvarives base function for BB base, key is a field name
Definition: EntitiesFieldData.cpp:869
MoFEM::EntitiesFieldData::EntData::getFTensor3Diff2N
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3Diff2N()
Get second derivatives of base functions for Hvec space.
Definition: EntitiesFieldData.hpp:779
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
HVEC0_2
@ HVEC0_2
Definition: definitions.h:198
MoFEM::EntitiesFieldData::EntData::getOrder
ApproximationOrder getOrder() const
get approximation order
Definition: EntitiesFieldData.hpp:1197
MoFEM::EntitiesFieldData::EntData::sPace
FieldSpace sPace
Entity space.
Definition: EntitiesFieldData.hpp:1068
MoFEM::constructor_derived_data
static void constructor_derived_data(DerivedEntitiesFieldData *derived_data, const boost::shared_ptr< EntitiesFieldData > &data_ptr)
Definition: EntitiesFieldData.cpp:86
NOBASE
@ NOBASE
Definition: definitions.h:59
MoFEM::EntitiesFieldData::EntData::getLocalIndices
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
Definition: EntitiesFieldData.hpp:1216
MoFEM::DerivedEntitiesFieldData::DerivedEntData::DerivedEntData
DerivedEntData(const boost::shared_ptr< EntitiesFieldData::EntData > &ent_data_ptr)
Definition: EntitiesFieldData.cpp:189
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::EntitiesFieldData::EntData::getBBNByOrderSharedPtr
virtual boost::shared_ptr< MatrixDouble > & getBBNByOrderSharedPtr(const size_t o)
get BB base by order
Definition: EntitiesFieldData.cpp:879
HVEC0_1
@ HVEC0_1
Definition: definitions.h:195
MoFEM::EntitiesFieldData::baseSwap
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Swap approximation base.
Definition: EntitiesFieldData.cpp:156
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1241
MoFEM::DerivedEntitiesFieldData::DerivedEntData::getBBDiffNSharedPtr
boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
Definition: EntitiesFieldData.cpp:928
MoFEM::EntitiesFieldData::EntData::getFTensor2N
auto getFTensor2N()
Get base functions for Hdiv space.
Definition: EntitiesFieldData.hpp:1521
MoFEM::EntitiesFieldData::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN()
Get derivatives of base functions for Hdiv space.
Definition: EntitiesFieldData.hpp:751
MoFEM::DerivedEntitiesFieldData
this class derive data form other data structure
Definition: EntitiesFieldData.hpp:109
MoFEM::EntitiesFieldData::EntData::getBBNMap
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBNMap()
get hash map of base function for BB base, key is a field name
Definition: EntitiesFieldData.cpp:864
MoFEM::EntitiesFieldData::EntData::getBBDiffNByOrderArray
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
Definition: EntitiesFieldData.cpp:902
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
HVEC0_0
@ HVEC0_0
Definition: definitions.h:192
HVEC1_1
@ HVEC1_1
Definition: definitions.h:196
MoFEM::EntitiesFieldData::EntData::getBBAlphaIndicesByOrderArray
virtual std::array< boost::shared_ptr< MatrixInt >, MaxBernsteinBezierOrder > & getBBAlphaIndicesByOrderArray()
Definition: EntitiesFieldData.cpp:890
MoFEM::EntitiesFieldData::EntitiesFieldData
EntitiesFieldData()
Definition: EntitiesFieldData.hpp:95
MoFEM::DerivedEntitiesFieldData::DerivedEntData::getBBAlphaIndicesSharedPtr
boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
Definition: EntitiesFieldData.cpp:907
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
MoFEM::DerivedEntitiesFieldData::DerivedEntData::getSense
int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
Definition: EntitiesFieldData.cpp:193
HVEC1
@ HVEC1
Definition: definitions.h:186
MoFEM::DerivedEntitiesFieldData::DerivedEntData
Derived ata on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:1127
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::EntitiesFieldData::EntData::getEntDataBitRefLevel
virtual std::vector< BitRefLevel > & getEntDataBitRefLevel()
Definition: EntitiesFieldData.cpp:944
FTensor::Tensor2::data
T data[Tensor_Dim0][Tensor_Dim1]
Definition: Tensor2_value.hpp:18
FTensor::Tensor1::data
T data[Tensor_Dim]
Definition: Tensor1_value.hpp:10
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
HVEC2_1
@ HVEC2_1
Definition: definitions.h:197
MoFEM::EntitiesFieldData::EntData::getBBNSharedPtr
virtual boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
Definition: EntitiesFieldData.cpp:828
MoFEM::EntitiesFieldData::EntData::getFTensor1N
auto getFTensor1N()
Get base functions for Hdiv space.
Definition: EntitiesFieldData.hpp:1511
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::EntitiesFieldData::EntData::MaxBernsteinBezierOrder
static constexpr size_t MaxBernsteinBezierOrder
Definition: EntitiesFieldData.hpp:1036
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
HVEC2_2
@ HVEC2_2
Definition: definitions.h:200
HVEC1_2
@ HVEC1_2
Definition: definitions.h:199
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::DerivedEntitiesFieldData::DerivedEntData::getNSharedPtr
boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives derivative)
Definition: EntitiesFieldData.cpp:198
MoFEM::EntitiesFieldData::EntData::BaseDerivatives
BaseDerivatives
Definition: EntitiesFieldData.hpp:129
MoFEM::EntitiesFieldData::EntData::bAse
FieldApproximationBase bAse
Field approximation base.
Definition: EntitiesFieldData.hpp:1069
MoFEM::EntitiesFieldData::EntData::iNdices
VectorInt iNdices
Global indices on entity.
Definition: EntitiesFieldData.hpp:1070
MoFEM::EntitiesFieldData::EntData::EntData
EntData(const bool allocate_base_matrices=true)
Definition: EntitiesFieldData.cpp:10
MoFEM::EntitiesFieldData::EntData::fieldEntities
VectorFieldEntities fieldEntities
Field entities.
Definition: EntitiesFieldData.hpp:1073
convert.type
type
Definition: convert.py:64
MoFEM::DerivedEntitiesFieldData::setElementType
MoFEMErrorCode setElementType(const EntityType type)
Definition: EntitiesFieldData.cpp:109
MoFEM::DerivedEntitiesFieldData::DerivedEntData::getDiffNSharedPtr
boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
Definition: EntitiesFieldData.cpp:216
MoFEM::EntitiesFieldData::EntData::getBBDiffNSharedPtr
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
Definition: EntitiesFieldData.cpp:845
MoFEM::DerivedEntitiesFieldData::dataPtr
const boost::shared_ptr< EntitiesFieldData > dataPtr
Definition: EntitiesFieldData.hpp:118
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
MoFEM::EntitiesFieldData::EntData::getBBNByOrderArray
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBNByOrderArray()
Definition: EntitiesFieldData.cpp:896
MoFEM::EntitiesFieldData::EntData::getDiffNSharedPtr
virtual boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
Definition: EntitiesFieldData.cpp:37
MoFEM::EntitiesFieldData::EntData::getBBAlphaIndicesByOrderSharedPtr
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesByOrderSharedPtr(const size_t o)
get ALpha indices for BB base by order
Definition: EntitiesFieldData.cpp:874
MoFEM::DerivedEntitiesFieldData::DerivedEntData::getBBNSharedPtr
boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
Definition: EntitiesFieldData.cpp:913
MoFEM::EntitiesFieldData::setElementType
virtual MoFEMErrorCode setElementType(const EntityType type)
Definition: EntitiesFieldData.cpp:79
MoFEM::EntitiesFieldData::operator<<
friend std::ostream & operator<<(std::ostream &os, const EntitiesFieldData &e)
Definition: EntitiesFieldData.cpp:270
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
MoFEM::EntitiesFieldData::bAse
std::bitset< LASTBASE > bAse
bases on element
Definition: EntitiesFieldData.hpp:45
MoFEM::constructor_data
static void constructor_data(EntitiesFieldData *data, const EntityType type)
Definition: EntitiesFieldData.cpp:42
MoFEM::EntitiesFieldData::EntData::getBBAlphaIndicesMap
virtual std::map< std::string, boost::shared_ptr< MatrixInt > > & getBBAlphaIndicesMap()
Definition: EntitiesFieldData.cpp:859
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN()
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:536
N
const int N
Definition: speed_test.cpp:3
MoFEM::EntitiesFieldData::EntData::getFTensor0FieldData
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
Definition: EntitiesFieldData.cpp:506
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
FTensor::Tensor0
Definition: Tensor0.hpp:16
MoFEM::EntitiesFieldData::EntData::localIndices
VectorInt localIndices
Local indices on entity.
Definition: EntitiesFieldData.hpp:1071
MoFEM::EntitiesFieldData::resetFieldDependentData
MoFEMErrorCode resetFieldDependentData()
Definition: EntitiesFieldData.cpp:127
MoFEM::EntitiesFieldData::EntData::getSense
virtual int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
Definition: EntitiesFieldData.hpp:1231
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::DerivedEntitiesFieldData::DerivedEntData::getEntDataBitRefLevel
std::vector< BitRefLevel > & getEntDataBitRefLevel()
Definition: EntitiesFieldData.cpp:949
MoFEM::EntitiesFieldData::EntData::getBBDiffNByOrderSharedPtr
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNByOrderSharedPtr(const size_t o)
get BB base derivative by order
Definition: EntitiesFieldData.cpp:884
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::EntitiesFieldData::EntData::LastDerivative
@ LastDerivative
Definition: EntitiesFieldData.hpp:135
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::DerivedEntitiesFieldData::DerivedEntData::baseSwap
MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Used by Bernstein base to keep temporally pointer.
Definition: EntitiesFieldData.cpp:168
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
HVEC2_0
@ HVEC2_0
Definition: definitions.h:194
MoFEM::EntitiesFieldData::EntData::baseFunctionsAndBaseDerivatives
std::array< std::array< boost::shared_ptr< MatrixDouble >, LASTBASE >, LastDerivative > baseFunctionsAndBaseDerivatives
Definition: EntitiesFieldData.hpp:1079
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:56
MoFEM::EntitiesFieldData::EntData::baseSwap
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Swap bases functions.
Definition: EntitiesFieldData.cpp:136
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::EntitiesFieldData::EntData::fieldData
VectorDouble fieldData
Field data on entity.
Definition: EntitiesFieldData.hpp:1074
MoFEM::DerivedEntitiesFieldData::DerivedEntitiesFieldData
DerivedEntitiesFieldData(const boost::shared_ptr< EntitiesFieldData > &data_ptr)
Definition: EntitiesFieldData.cpp:103
HVEC2
@ HVEC2
Definition: definitions.h:186
MoFEM::EntitiesFieldData::EntData::dOfs
VectorDofs dOfs
DoFs on entity.
Definition: EntitiesFieldData.hpp:1072
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
HVEC1_0
@ HVEC1_0
Definition: definitions.h:193
HVEC0
@ HVEC0
Definition: definitions.h:186
MoFEM::EntitiesFieldData::EntData::getNSharedPtr
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives direvatie)
Definition: EntitiesFieldData.cpp:27