v0.10.0
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  // FIXME: precision should not be set here
280  os << "fieldData: " << std::fixed << std::setprecision(2) << 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 
static Index< 'o', 3 > o
T data[Tensor_Dim]
T data[Tensor_Dim0][Tensor_Dim1]
FieldApproximationBase
approximation base
Definition: definitions.h:150
@ LASTBASE
Definition: definitions.h:161
@ NOBASE
Definition: definitions.h:151
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
@ NOSPACE
Definition: definitions.h:175
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ HVEC0
Definition: definitions.h:255
@ HVEC1
Definition: definitions.h:255
@ HVEC2
Definition: definitions.h:255
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
@ HVEC1_1
Definition: definitions.h:265
@ HVEC0_1
Definition: definitions.h:264
@ HVEC1_0
Definition: definitions.h:262
@ HVEC2_1
Definition: definitions.h:266
@ HVEC1_2
Definition: definitions.h:268
@ HVEC2_2
Definition: definitions.h:269
@ HVEC2_0
Definition: definitions.h:263
@ HVEC0_2
Definition: definitions.h:267
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
static void constructor_derived_data(DerivedDataForcesAndSourcesCore *derived_data, const boost::shared_ptr< DataForcesAndSourcesCore > &data_ptr)
static void constructor_data(DataForcesAndSourcesCore *data, const EntityType type)
DataForcesAndSourcesCore::EntData EntData
const int N
Definition: speed_test.cpp:3
Data on single entity (This is passed as argument to DataOperator::doWork)
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBNMap()
virtual boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
std::array< boost::shared_ptr< MatrixDouble >, LASTBASE > diffN
Derivatives of base functions.
virtual boost::shared_ptr< MatrixDouble > & getBBNByOrderSharedPtr(const size_t o)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
FieldApproximationBase bAse
Field approximation base.
VectorInt iNdices
Global indices on entity.
virtual int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
virtual boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN()
Get derivatives of base functions.
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
std::array< boost::shared_ptr< MatrixDouble >, LASTBASE > N
Base functions.
auto getFTensor2N()
Get base functions for Hdiv space.
VectorInt localIndices
Local indices on entity.
virtual std::array< boost::shared_ptr< MatrixInt >, MaxBernsteinBezierOrder > & getBBAlphaIndicesByOrderArray()
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
auto getFTensor1N()
Get base functions for Hdiv space.
const VectorDouble & getFieldData() const
get dofs values
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesByOrderSharedPtr(const size_t o)
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN()
Get derivatives of base functions for Hdiv space.
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBNByOrderArray()
VectorDouble fieldData
Field data on entity.
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
EntData(const bool allocate_base_matrices=true)
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
virtual std::map< std::string, boost::shared_ptr< MatrixInt > > & getBBAlphaIndicesMap()
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNByOrderSharedPtr(const size_t o)
ApproximationOrder getOrder() const
get approximation order
static constexpr size_t MaxBernsteinBezierOrder
data structure for finite element entity
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Swap approximation base.
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MoFEMErrorCode resetFieldDependentData()
std::bitset< LASTBASE > bAse
bases on element
friend std::ostream & operator<<(std::ostream &os, const DataForcesAndSourcesCore &e)
virtual MoFEMErrorCode setElementType(const EntityType type)
Derived ata on single entity (This is passed as argument to DataOperator::doWork)
boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base)
int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
DerivedEntData(const boost::shared_ptr< DataForcesAndSourcesCore::EntData > &ent_data_ptr)
this class derive data form other data structure
const boost::shared_ptr< DataForcesAndSourcesCore > dataPtr
MoFEMErrorCode setElementType(const EntityType type)
DerivedDataForcesAndSourcesCore(const boost::shared_ptr< DataForcesAndSourcesCore > &data_ptr)