v0.9.1
DataStructures.cpp
Go to the documentation of this file.
1 /** \file DataStructures.cpp
2 \brief Implementation for Data Structures in Forces and Sources
3 
4 */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 namespace MoFEM {
21 
22 DataForcesAndSourcesCore::EntData::EntData(const bool allocate_base_matrices)
23  : sEnse(0), oRder(0), bAse(NOBASE) {
24  if (allocate_base_matrices)
25  for (int b = 0; b != LASTBASE; ++b) {
26  N[b].reset(new MatrixDouble());
27  diffN[b].reset(new MatrixDouble());
28  }
29 }
30 
31 int DataForcesAndSourcesCore::EntData::getSense() const { return sEnse; }
32 
33 boost::shared_ptr<MatrixDouble> &
35  const FieldApproximationBase base) {
36  return N[base];
37 }
38 
39 const boost::shared_ptr<MatrixDouble> &
41  const FieldApproximationBase base) const {
42  return N[base];
43 }
44 
45 boost::shared_ptr<MatrixDouble> &
47  const FieldApproximationBase base) {
48  return diffN[base];
49 }
50 
51 const boost::shared_ptr<MatrixDouble> &
53  const FieldApproximationBase base) const {
54  return diffN[base];
55 }
56 
58  const EntityType type) {
59 
61 
62  data->dataOnEntities[MBENTITYSET].push_back(new EntData());
63 
64  switch (type) {
65  case MBENTITYSET:
66  break;
67  case MBTET:
68  data->dataOnEntities[MBVERTEX].push_back(new EntData());
69  for (int ee = 0; ee != 6; ++ee) {
70  data->dataOnEntities[MBEDGE].push_back(new EntData());
71  }
72  for (int ff = 0; ff != 4; ++ff) {
73  data->dataOnEntities[MBTRI].push_back(new EntData());
74  }
75  data->dataOnEntities[MBTET].push_back(new EntData());
76  break;
77  case MBTRI:
78  data->dataOnEntities[MBVERTEX].push_back(new EntData());
79  for (int ee = 0; ee != 3; ++ee) {
80  data->dataOnEntities[MBEDGE].push_back(new EntData());
81  }
82  data->dataOnEntities[MBTRI].push_back(new EntData());
83  break;
84  case MBQUAD:
85  data->dataOnEntities[MBVERTEX].push_back(new EntData());
86  for (int ee = 0; ee != 4; ++ee) {
87  data->dataOnEntities[MBEDGE].push_back(new EntData());
88  }
89  data->dataOnEntities[MBQUAD].push_back(new EntData());
90  break;
91  case MBEDGE:
92  data->dataOnEntities[MBVERTEX].push_back(new EntData());
93  data->dataOnEntities[MBEDGE].push_back(new EntData());
94  break;
95  case MBVERTEX:
96  data->dataOnEntities[MBVERTEX].push_back(new EntData());
97  break;
98  case MBPRISM:
99  data->dataOnEntities[MBVERTEX].push_back(new EntData());
100  for (int ee = 0; ee != 9; ++ee) {
101  data->dataOnEntities[MBEDGE].push_back(new EntData());
102  }
103  for (int qq = 0; qq != 5; ++qq) {
104  data->dataOnEntities[MBQUAD].push_back(new EntData());
105  }
106  for (int ff = 0; ff != 5; ++ff) {
107  data->dataOnEntities[MBTRI].push_back(new EntData());
108  }
109  data->dataOnEntities[MBPRISM].push_back(new EntData());
110  break;
111  default:
113  }
114 }
115 
117  constructor_data(this, type);
118 }
119 
122  for (auto &data : dataOnEntities)
123  data.clear();
124  constructor_data(this, type);
126 }
127 
129  DerivedDataForcesAndSourcesCore *derived_data,
130  const boost::shared_ptr<DataForcesAndSourcesCore> &data_ptr) {
131 
134 
135  for (int tt = MBVERTEX; tt != MBMAXTYPE; ++tt) {
136  auto &ent_data = data_ptr->dataOnEntities[tt];
137  auto &derived_ent_data = derived_data->dataOnEntities[tt];
138  for (auto &e : ent_data) {
139  boost::shared_ptr<EntData> ent_data_ptr(data_ptr, &e);
140  derived_ent_data.push_back(new DerivedEntData(ent_data_ptr));
141  }
142  }
143 }
144 
146  const boost::shared_ptr<DataForcesAndSourcesCore> &data_ptr)
147  : DataForcesAndSourcesCore(), dataPtr(data_ptr) {
149 }
150 
154  for (EntityType tt = MBVERTEX; tt != MBMAXTYPE; ++tt)
155  dataOnEntities[tt].clear();
158 }
159 
162  sPace = NOSPACE;
163  bAse = NOBASE;
164  iNdices.resize(0, false);
165  localIndices.resize(0, false);
166  dOfs.resize(0, false);
167  fieldData.resize(0, false);
169 }
170 
173  for (EntityType t = MBVERTEX; t != MBMAXTYPE; t++)
174  for (auto &e : dataOnEntities[t])
175  CHKERR e.resetFieldDependentData();
177 }
178 
180 DataForcesAndSourcesCore::EntData::baseSwap(const std::string &field_name,
181  const FieldApproximationBase base) {
183  auto make_swap = [](boost::shared_ptr<MatrixDouble> &ptr,
184  boost::shared_ptr<MatrixDouble> &ptrBB,
185  boost::shared_ptr<MatrixDouble> &swap_ptr) {
186  if (swap_ptr) {
187  ptr = swap_ptr;
188  swap_ptr.reset();
189  } else {
190  swap_ptr = ptr;
191  ptr = ptrBB;
192  }
193  };
194  make_swap(getNSharedPtr(base), getBBNSharedPtr(field_name), swapBaseNPtr);
195  make_swap(getDiffNSharedPtr(base), getBBDiffNSharedPtr(field_name),
196  swapBaseDiffNPtr);
198 }
199 
201  const FieldApproximationBase base) {
203  for (int tt = MBVERTEX; tt != MBMAXTYPE; ++tt) {
204  auto &ent_data = dataOnEntities[tt];
205  for (auto &side_data : ent_data)
206  CHKERR side_data.baseSwap(field_name, base);
207  }
209 }
210 
212  const std::string &field_name, const FieldApproximationBase base) {
214  auto make_swap = [](boost::shared_ptr<MatrixDouble> &ptr,
215  boost::shared_ptr<MatrixDouble> &ptrBB,
216  boost::shared_ptr<MatrixDouble> &swap_ptr) {
217  if (swap_ptr) {
218  ptr = swap_ptr;
219  swap_ptr.reset();
220  } else {
221  swap_ptr = ptr;
222  ptr = ptrBB;
223  }
224  };
225  make_swap(getDerivedNSharedPtr(base), getBBNSharedPtr(field_name),
226  swapBaseNPtr);
227  make_swap(getDerivedDiffNSharedPtr(base), getBBDiffNSharedPtr(field_name),
228  swapBaseDiffNPtr);
230 }
231 
233  const boost::shared_ptr<DataForcesAndSourcesCore::EntData> &ent_data_ptr)
234  : DataForcesAndSourcesCore::EntData(false), entDataPtr(ent_data_ptr) {}
235 
237  return entDataPtr->getSense();
238 }
239 
240 boost::shared_ptr<MatrixDouble> &
242  const FieldApproximationBase base) {
243  if (N[base])
244  return N[base];
245  else
246  return entDataPtr->getNSharedPtr(base);
247 }
248 boost::shared_ptr<MatrixDouble> &
250  const FieldApproximationBase base) {
251  if (diffN[base])
252  return diffN[base];
253  else
254  return entDataPtr->getDiffNSharedPtr(base);
255 }
256 const boost::shared_ptr<MatrixDouble> &
258  const FieldApproximationBase base) const {
259  if (N[base])
260  return N[base];
261  else
262  return entDataPtr->getNSharedPtr(base);
263 }
264 const boost::shared_ptr<MatrixDouble> &
266  const FieldApproximationBase base) const {
267  if (diffN[base])
268  return diffN[base];
269  else
270  return entDataPtr->getDiffNSharedPtr(base);
271 }
272 
273 std::ostream &operator<<(std::ostream &os,
275  os << "sEnse: " << e.getSense() << std::endl
276  << "oRder: " << e.getOrder() << std::endl
277  << "global indices: " << e.getIndices() << std::endl
278  << "local indices: " << e.getLocalIndices() << std::endl;
279  os.precision(2);
280  os << "fieldData: " << std::fixed << e.getFieldData() << std::endl;
281  MatrixDouble base = const_cast<DataForcesAndSourcesCore::EntData &>(e).getN();
282  MatrixDouble diff_base =
283  const_cast<DataForcesAndSourcesCore::EntData &>(e).getDiffN();
284  const double eps = 1e-6;
285  for (unsigned int ii = 0; ii != base.size1(); ii++) {
286  for (unsigned int jj = 0; jj != base.size2(); jj++) {
287  if (fabs(base(ii, jj)) < eps)
288  base(ii, jj) = 0;
289  }
290  }
291  for (unsigned int ii = 0; ii != diff_base.size1(); ii++) {
292  for (unsigned int jj = 0; jj != diff_base.size2(); jj++) {
293  if (fabs(diff_base(ii, jj)) < eps)
294  diff_base(ii, jj) = 0;
295  }
296  }
297  os << "N: " << std::fixed << base << std::endl
298  << "diffN: " << std::fixed << diff_base;
299  return os;
300 }
301 
302 std::ostream &operator<<(std::ostream &os, const DataForcesAndSourcesCore &e) {
303  for (unsigned int nn = 0; nn < e.dataOnEntities[MBVERTEX].size(); nn++) {
304  os << "dataOnEntities[MBVERTEX][" << nn << "]" << std::endl
305  << e.dataOnEntities[MBVERTEX][nn] << std::endl;
306  }
307  for (unsigned int ee = 0; ee < e.dataOnEntities[MBEDGE].size(); ee++) {
308  os << "dataOnEntities[MBEDGE][" << ee << "]" << std::endl
309  << e.dataOnEntities[MBEDGE][ee] << std::endl;
310  }
311  for (unsigned int ff = 0; ff < e.dataOnEntities[MBTRI].size(); ff++) {
312  os << "dataOnEntities[MBTRI][" << ff << "] " << std::endl
313  << e.dataOnEntities[MBTRI][ff] << std::endl;
314  }
315  for (unsigned int vv = 0; vv < e.dataOnEntities[MBTET].size(); vv++) {
316  os << "dataOnEntities[MBTET][" << vv << "]" << std::endl
317  << e.dataOnEntities[MBTET][vv] << std::endl;
318  }
319  return os;
320 }
321 
322 /** \name Specializations for H1/L2 */
323 
324 /**@{*/
325 
326 template <>
328 DataForcesAndSourcesCore::EntData::getFTensor1FieldData<3>() {
329  if (dOfs[0]->getNbOfCoeffs() != 3) {
330  std::stringstream s;
331  s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
332  s << " but you ask for tensor rank 1 dimension 3";
333  THROW_MESSAGE(s.str());
334  }
335  double *ptr = &*fieldData.data().begin();
336  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
337  &ptr[2]);
338 }
339 
340 template <>
342 DataForcesAndSourcesCore::EntData::getFTensor1FieldData<2>() {
343  if (dOfs[0]->getNbOfCoeffs() != 2) {
344  std::stringstream s;
345  s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
346  s << " but you ask for tensor rank 1 dimension 3";
347  THROW_MESSAGE(s.str());
348  }
349  double *ptr = &*fieldData.data().begin();
350  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(ptr, &ptr[1]);
351 }
352 
353 template <>
355 DataForcesAndSourcesCore::EntData::getFTensor2FieldData<3, 3>() {
356  if (dOfs[0]->getNbOfCoeffs() != 9) {
357  std::stringstream s;
358  s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
359  s << " but you ask for tensor rank 2 dimensions 3 by 3 so 9 coefficients "
360  "is expected";
361  THROW_MESSAGE(s.str());
362  }
363  double *ptr = &*fieldData.data().begin();
365  ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
366  &ptr[8]);
367 }
368 
369 template <>
371 DataForcesAndSourcesCore::EntData::getFTensor2SymmetricFieldData<3>() {
372  if (dOfs[0]->getNbOfCoeffs() != 6) {
373  std::stringstream s;
374  s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
375  s << " but you ask for symmetric tensor rank 2 dimensions 3 by 3 so 6 "
376  "coefficients "
377  "is expected";
378  THROW_MESSAGE(s.str());
379  }
380  double *ptr = &*fieldData.data().begin();
382  ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
383 }
384 
385 template <>
387 DataForcesAndSourcesCore::EntData::getFTensor2SymmetricFieldData<2>() {
388  if (dOfs[0]->getNbOfCoeffs() != 3) {
389  std::stringstream s;
390  s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
391  s << " but you ask for symmetric tensor rank 2 dimensions 2 by 2 so 3 "
392  "coefficients "
393  "is expected";
394  THROW_MESSAGE(s.str());
395  }
396  double *ptr = &*fieldData.data().begin();
398  ptr, &ptr[1], &ptr[2]);
399 }
400 
403  if (dOfs[0]->getNbOfCoeffs() != 1) {
404  std::stringstream s;
405  s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
406  s << " but expected scalar field, tensor of rank 0";
407  THROW_MESSAGE(s.str());
408  }
410  &*fieldData.data().begin());
411 }
412 
413 template <int Tensor_Dim>
416  const FieldApproximationBase base) {
417  std::stringstream s;
418  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
419  THROW_MESSAGE(s.str());
421 }
422 
423 template <int Tensor_Dim>
426  return getFTensor1DiffN<Tensor_Dim>(bAse);
427 }
428 
429 template <int Tensor_Dim>
432  const FieldApproximationBase base, const int bb) {
433  std::stringstream s;
434  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
435  THROW_MESSAGE(s.str());
437 }
438 
439 template <int Tensor_Dim>
442  return getFTensor1DiffN<Tensor_Dim>(bAse, bb);
443 }
444 
445 /**
446  * \brief Get spatial derivative of base function tensor for dimension 3d
447  */
448 template <>
450 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
451  const FieldApproximationBase base) {
452  double *ptr = &*getDiffN(base).data().begin();
453  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
454 }
455 
456 /**
457  * \brief Get spatial derivative of base function tensor for dimension 3d
458  */
459 template <>
461 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>() {
462  return getFTensor1DiffN<3>(bAse);
463 }
464 
465 /**
466  * \brief Get spatial derivative of base function tensor for dimension 3d
467  */
468 template <>
470 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
471  const FieldApproximationBase base, const int bb) {
472  double *ptr = &*getDiffN(base).data().begin();
474  &ptr[3 * bb], &ptr[3 * bb + 1], &ptr[3 * bb + 2], getDiffN(base).size2());
475 }
476 
477 /**
478  * \brief Get spatial derivative of base function tensor for dimension 3d
479  */
480 template <>
482 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(const int bb) {
483  return getFTensor1DiffN<3>(bAse, bb);
484 }
485 
486 /**
487  * \brief Get spatial derivative of base function tensor for dimension 2d
488  */
489 template <>
491 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
492  const FieldApproximationBase base) {
493  double *ptr = &*getDiffN(base).data().begin();
494  return FTensor::Tensor1<double *, 2>(ptr, &ptr[1], 2);
495 }
496 
497 /**
498  * \brief Get spatial derivative of base function tensor for dimension 2d
499  */
500 template <>
502 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>() {
503  return getFTensor1DiffN<2>(bAse);
504 }
505 
506 /**
507  * \brief Get spatial derivative of base function tensor for dimension 3d
508  */
509 template <>
511 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
512  const FieldApproximationBase base, const int bb) {
513  double *ptr = &*getDiffN(base).data().begin();
514  return FTensor::Tensor1<double *, 2>(&ptr[2 * bb], &ptr[2 * bb + 1],
515  getDiffN(base).size1());
516 }
517 
518 /**
519  * \brief Get spatial derivative of base function tensor for dimension 3d
520  */
521 template <>
523 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(const int bb) {
524  return getFTensor1DiffN<2>(bAse, bb);
525 }
526 
527 template <int Tensor_Dim>
530  const FieldApproximationBase base, const int gg, const int bb) {
531  std::stringstream s;
532  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
533  THROW_MESSAGE(s.str());
535 }
536 
537 /**
538  * \brief Get spatial derivative of base function tensor for dimension 3d
539  */
540 template <>
542 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
543  const FieldApproximationBase base, const int gg, const int bb) {
544  double *ptr = &getDiffN(base)(gg, 3 * bb);
545  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
546 }
547 
548 /**
549  * \brief Get spatial derivative of base function tensor for dimension 3d
550  */
551 template <>
553 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(const int gg,
554  const int bb) {
555  return getFTensor1DiffN<3>(bAse, gg, bb);
556 }
557 
558 /**
559  * \brief Get spatial derivative of base function tensor for dimension 2d
560  */
561 template <>
563 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
564  const FieldApproximationBase base, const int gg, const int bb) {
565  double *ptr = &getDiffN(base)(gg, 2 * bb);
566  return FTensor::Tensor1<double *, 2>(ptr, &ptr[1], 2);
567 }
568 
569 /**
570  * \brief Get spatial derivative of base function tensor for dimension 2d
571  */
572 template <>
574 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(const int gg,
575  const int bb) {
576  return getFTensor1DiffN<2>(bAse, gg, bb);
577 }
578 
579 /**@}*/
580 
581 /** \name Specializations for HDiv/HCrul */
582 
583 /**@{*/
584 
585 template <int Tensor_Dim>
588  std::stringstream s;
589  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
590  THROW_MESSAGE(s.str());
592 }
593 
594 template <int Tensor_Dim>
597  const int gg, const int bb) {
598  std::stringstream s;
599  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
600  THROW_MESSAGE(s.str());
602 }
603 
604 template <int Tensor_Dim0, int Tensor_Dim1>
606  Tensor_Dim0, Tensor_Dim1>
608  FieldApproximationBase base) {
609  std::stringstream s;
610  s << "Template for tensor dimension " << Tensor_Dim0 << "x" << Tensor_Dim1
611  << " not implemented";
612  THROW_MESSAGE(s.str());
614 }
615 
616 template <int Tensor_Dim0, int Tensor_Dim1>
618  Tensor_Dim0, Tensor_Dim1>
620  const int gg,
621  const int bb) {
622  std::stringstream s;
623  s << "Template for tensor dimension " << Tensor_Dim0 << "x" << Tensor_Dim1
624  << " not implemented";
625  THROW_MESSAGE(s.str());
627 }
628 
629 template <>
631 DataForcesAndSourcesCore::EntData::getFTensor1N<3>(
632  FieldApproximationBase base) {
633  double *t_n_ptr = &*getN(base).data().begin();
634  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(t_n_ptr, // HVEC0
635  &t_n_ptr[HVEC1],
636  &t_n_ptr[HVEC2]);
637 }
638 
639 template <>
641 DataForcesAndSourcesCore::EntData::getFTensor1N<3>(FieldApproximationBase base,
642  const int gg, const int bb) {
643  double *t_n_ptr = &getN(base)(gg, 3 * bb);
644  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(t_n_ptr, // HVEC0
645  &t_n_ptr[HVEC1],
646  &t_n_ptr[HVEC2]);
647 }
648 
649 template <>
651 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 3>(
652  FieldApproximationBase base) {
653  double *t_diff_n_ptr = &*getDiffN(base).data().begin();
655  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC0_2],
656  &t_diff_n_ptr[HVEC1_0], &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC1_2],
657  &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1], &t_diff_n_ptr[HVEC2_2]);
658 }
659 
660 template <>
662 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 3>(
663  FieldApproximationBase base, const int gg, const int bb) {
664  double *t_diff_n_ptr = &getDiffN(base)(gg, 9 * bb);
666  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC0_2],
667  &t_diff_n_ptr[HVEC1_0], &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC1_2],
668  &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1], &t_diff_n_ptr[HVEC2_2]);
669 }
670 
671 template <>
673 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 2>(
674  FieldApproximationBase base) {
675  double *t_diff_n_ptr = &*getDiffN(base).data().begin();
677  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC1_0],
678  &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1]);
679 }
680 
681 template <>
683 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 2>(
684  FieldApproximationBase base, const int gg, const int bb) {
685  double *t_diff_n_ptr = &getDiffN(base)(gg, 6 * bb);
687  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC1_0],
688  &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1]);
689 }
690 
691 template <int Tensor_Dim0, int Tensor_Dim1>
693  Tensor_Dim0, Tensor_Dim1>
695  std::stringstream s;
696  s << "Template for tensor dimension " << Tensor_Dim0 << ", " << Tensor_Dim1
697  << " not implemented";
698  THROW_MESSAGE(s.str());
700  Tensor_Dim0, Tensor_Dim1>();
701 }
702 
703 template <>
705 DataForcesAndSourcesCore::EntData::getFTensor2N<3, 3>(
706  FieldApproximationBase base) {
707  double *t_n_ptr = &*(getN(base).data().begin());
709 
710  &t_n_ptr[HVEC0], &t_n_ptr[HVEC1], &t_n_ptr[HVEC2],
711 
712  &t_n_ptr[3 + HVEC0], &t_n_ptr[3 + HVEC1], &t_n_ptr[3 + HVEC2],
713 
714  &t_n_ptr[6 + HVEC0], &t_n_ptr[6 + HVEC1], &t_n_ptr[6 + HVEC2]
715 
716  );
717 }
718 
719 template <int Tensor_Dim0, int Tensor_Dim1>
721  Tensor_Dim0, Tensor_Dim1>
723  const int gg, const int bb) {
724  std::stringstream s;
725  s << "Template for tensor dimension " << Tensor_Dim0 << ", " << Tensor_Dim1
726  << " not implemented";
727  THROW_MESSAGE(s.str());
729  Tensor_Dim0, Tensor_Dim1>();
730 }
731 
732 template <>
734 DataForcesAndSourcesCore::EntData::getFTensor2N<3, 3>(
735  FieldApproximationBase base, const int gg, const int bb) {
736  double *t_n_ptr = &getN(base)(gg, 9 * bb);
738 
739  &t_n_ptr[HVEC0], &t_n_ptr[HVEC1], &t_n_ptr[HVEC2],
740 
741  &t_n_ptr[3 + HVEC0], &t_n_ptr[3 + HVEC1], &t_n_ptr[3 + HVEC2],
742 
743  &t_n_ptr[6 + HVEC0], &t_n_ptr[6 + HVEC1], &t_n_ptr[6 + HVEC2]
744 
745  );
746 }
747 
748 /**@}*/
749 
750 /** \name Bernstein-Bezier base only functions */
751 
752 /**@{*/
753 
754 boost::shared_ptr<MatrixInt> &
756  const std::string &field_name) {
757  return bbAlphaInduces[field_name];
758 }
759 
760 boost::shared_ptr<MatrixDouble> &
762  const std::string &field_name) {
763  return bbN[field_name];
764 }
765 
766 /**
767  * Get shared pointer to BB base base functions
768  */
769 const boost::shared_ptr<MatrixDouble> &
771  const std::string &field_name) const {
772  return bbN.at(field_name);
773 }
774 
775 /**
776  * Get shared pointer to BB derivatives of base base functions
777  */
778 boost::shared_ptr<MatrixDouble> &
780  const std::string &field_name) {
781  return bbDiffN[field_name];
782 }
783 
784 /**
785  * Get shared pointer to derivatives of BB base base functions
786  */
787 const boost::shared_ptr<MatrixDouble> &
789  const std::string &field_name) const {
790  return bbDiffN.at(field_name);
791 }
792 
793 std::map<std::string, boost::shared_ptr<MatrixInt>> &
795  return bbAlphaInduces;
796 }
797 
798 std::map<std::string, boost::shared_ptr<MatrixDouble>> &
800  return bbN;
801 }
802 
803 std::map<std::string, boost::shared_ptr<MatrixDouble>> &
805  return bbDiffN;
806 }
807 
808 boost::shared_ptr<MatrixInt> &
810  const size_t o) {
811  return bbAlphaInducesByOrder[o];
812 }
813 
814 boost::shared_ptr<MatrixDouble> &
816  return bbNByOrder[o];
817 }
818 
819 const boost::shared_ptr<MatrixDouble> &
821  return bbNByOrder[o];
822 }
823 
824 boost::shared_ptr<MatrixDouble> &
826  return bbDiffNByOrder[o];
827 }
828 
829 const boost::shared_ptr<MatrixDouble> &
831  const size_t o) const {
832  return bbDiffNByOrder[o];
833 }
834 
835 std::array<boost::shared_ptr<MatrixInt>,
838  return bbAlphaInducesByOrder;
839 }
840 
841 std::array<boost::shared_ptr<MatrixDouble>,
844  return bbNByOrder;
845 }
846 
847 std::array<boost::shared_ptr<MatrixDouble>,
850  return bbDiffNByOrder;
851 }
852 
853 boost::shared_ptr<MatrixInt> &
855  const std::string &field_name) {
856  return entDataPtr->getBBAlphaIndicesSharedPtr(field_name);
857 }
858 
859 boost::shared_ptr<MatrixDouble> &
861  const std::string &field_name) {
862  return entDataPtr->getBBNSharedPtr(field_name);
863 }
864 
865 const boost::shared_ptr<MatrixDouble> &
867  const std::string &field_name) const {
868  return entDataPtr->getBBNSharedPtr(field_name);
869 }
870 
871 /**
872  * Get shared pointer to BB derivatives of base base functions
873  */
874 boost::shared_ptr<MatrixDouble> &
876  const std::string &field_name) {
877  return entDataPtr->getBBDiffNSharedPtr(field_name);
878 }
879 
880 /**
881  * Get shared pointer to derivatives of BB base base functions
882  */
883 const boost::shared_ptr<MatrixDouble> &
885  const std::string &field_name) const {
886  return entDataPtr->getBBDiffNSharedPtr(field_name);
887 }
888 
889 /**@}*/
890 
891 
892 } // namespace MoFEM
893 
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBNMap()
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
DerivedDataForcesAndSourcesCore(const boost::shared_ptr< DataForcesAndSourcesCore > &data_ptr)
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
static void constructor_derived_data(DerivedDataForcesAndSourcesCore *derived_data, const boost::shared_ptr< DataForcesAndSourcesCore > &data_ptr)
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
std::array< boost::shared_ptr< MatrixDouble >, LASTBASE > N
Base functions.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
VectorInt iNdices
Global indices on entity.
static void constructor_data(DataForcesAndSourcesCore *data, const EntityType type)
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN()
Get derivatives of base functions for Hdiv space.
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:625
MoFEMErrorCode setElementType(const EntityType type)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
virtual MoFEMErrorCode setElementType(const EntityType type)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
T data[Tensor_Dim]
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
virtual boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBNByOrderArray()
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesByOrderSharedPtr(const size_t o)
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
VectorDouble fieldData
Field data on entity.
FieldApproximationBase
approximation base
Definition: definitions.h:149
const boost::shared_ptr< DataForcesAndSourcesCore > dataPtr
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
virtual int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
virtual boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
static constexpr size_t MaxBernsteinBezierOrder
boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
virtual boost::shared_ptr< MatrixDouble > & getBBNByOrderSharedPtr(const size_t o)
ApproximationOrder getOrder() const
get approximation order
#define CHKERR
Inline error check.
Definition: definitions.h:601
Derived ata on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base)
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Swap approximation base.
FieldApproximationBase bAse
Field approximation base.
EntData(const bool allocate_base_matrices=true)
T data[Tensor_Dim0][Tensor_Dim1]
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN()
Get derivatives of base functions.
friend std::ostream & operator<<(std::ostream &os, const DataForcesAndSourcesCore &e)
DerivedEntData(const boost::shared_ptr< DataForcesAndSourcesCore::EntData > &ent_data_ptr)
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
Data on single entity (This is passed as argument to DataOperator::doWork)
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNByOrderSharedPtr(const size_t o)
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base)
virtual std::map< std::string, boost::shared_ptr< MatrixInt > > & getBBAlphaIndicesMap()
auto getFTensor2N()
Get base functions for Hdiv space.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
static const double eps
MoFEMErrorCode resetFieldDependentData()
const VectorDouble & getFieldData() const
get dofs values
const int N
Definition: speed_test.cpp:3
std::bitset< LASTBASE > bAse
bases on element
virtual std::array< boost::shared_ptr< MatrixInt >, MaxBernsteinBezierOrder > & getBBAlphaIndicesByOrderArray()
VectorInt localIndices
Local indices on entity.
auto getFTensor1N()
Get base functions for Hdiv space.
this class derive data form other data structure
std::array< boost::shared_ptr< MatrixDouble >, LASTBASE > diffN
Derivatives of base functions.