v0.9.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  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 
387  if (dOfs[0]->getNbOfCoeffs() != 1) {
388  std::stringstream s;
389  s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
390  s << " but expected scalar field, tensor of rank 0";
391  THROW_MESSAGE(s.str());
392  }
394  &*fieldData.data().begin());
395 }
396 
397 template <int Tensor_Dim>
400  const FieldApproximationBase base) {
401  std::stringstream s;
402  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
403  THROW_MESSAGE(s.str());
405 }
406 
407 template <int Tensor_Dim>
410  return getFTensor1DiffN<Tensor_Dim>(bAse);
411 }
412 
413 template <int Tensor_Dim>
416  const FieldApproximationBase base, const int bb) {
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, bb);
427 }
428 
429 /**
430  * \brief Get spatial derivative of base function tensor for dimension 3d
431  */
432 template <>
434 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
435  const FieldApproximationBase base) {
436  double *ptr = &*getDiffN(base).data().begin();
437  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
438 }
439 
440 /**
441  * \brief Get spatial derivative of base function tensor for dimension 3d
442  */
443 template <>
445 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>() {
446  return getFTensor1DiffN<3>(bAse);
447 }
448 
449 /**
450  * \brief Get spatial derivative of base function tensor for dimension 3d
451  */
452 template <>
454 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
455  const FieldApproximationBase base, const int bb) {
456  double *ptr = &*getDiffN(base).data().begin();
458  &ptr[3 * bb], &ptr[3 * bb + 1], &ptr[3 * bb + 2], getDiffN(base).size2());
459 }
460 
461 /**
462  * \brief Get spatial derivative of base function tensor for dimension 3d
463  */
464 template <>
466 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(const int bb) {
467  return getFTensor1DiffN<3>(bAse, bb);
468 }
469 
470 /**
471  * \brief Get spatial derivative of base function tensor for dimension 2d
472  */
473 template <>
475 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
476  const FieldApproximationBase base) {
477  double *ptr = &*getDiffN(base).data().begin();
478  return FTensor::Tensor1<double *, 2>(ptr, &ptr[1], 2);
479 }
480 
481 /**
482  * \brief Get spatial derivative of base function tensor for dimension 2d
483  */
484 template <>
486 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>() {
487  return getFTensor1DiffN<2>(bAse);
488 }
489 
490 /**
491  * \brief Get spatial derivative of base function tensor for dimension 3d
492  */
493 template <>
495 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
496  const FieldApproximationBase base, const int bb) {
497  double *ptr = &*getDiffN(base).data().begin();
498  return FTensor::Tensor1<double *, 2>(&ptr[2 * bb], &ptr[2 * bb + 1],
499  getDiffN(base).size1());
500 }
501 
502 /**
503  * \brief Get spatial derivative of base function tensor for dimension 3d
504  */
505 template <>
507 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(const int bb) {
508  return getFTensor1DiffN<2>(bAse, bb);
509 }
510 
511 template <int Tensor_Dim>
514  const FieldApproximationBase base, const int gg, const int bb) {
515  std::stringstream s;
516  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
517  THROW_MESSAGE(s.str());
519 }
520 
521 /**
522  * \brief Get spatial derivative of base function tensor for dimension 3d
523  */
524 template <>
526 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
527  const FieldApproximationBase base, const int gg, const int bb) {
528  double *ptr = &getDiffN(base)(gg, 3 * bb);
529  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
530 }
531 
532 /**
533  * \brief Get spatial derivative of base function tensor for dimension 3d
534  */
535 template <>
537 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(const int gg,
538  const int bb) {
539  return getFTensor1DiffN<3>(bAse, gg, bb);
540 }
541 
542 /**
543  * \brief Get spatial derivative of base function tensor for dimension 2d
544  */
545 template <>
547 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
548  const FieldApproximationBase base, const int gg, const int bb) {
549  double *ptr = &getDiffN(base)(gg, 2 * bb);
550  return FTensor::Tensor1<double *, 2>(ptr, &ptr[1], 2);
551 }
552 
553 /**
554  * \brief Get spatial derivative of base function tensor for dimension 2d
555  */
556 template <>
558 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(const int gg,
559  const int bb) {
560  return getFTensor1DiffN<2>(bAse, gg, bb);
561 }
562 
563 /**@}*/
564 
565 /** \name Specializations for HDiv/HCrul */
566 
567 /**@{*/
568 
569 template <int Tensor_Dim>
572  std::stringstream s;
573  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
574  THROW_MESSAGE(s.str());
576 }
577 
578 template <int Tensor_Dim>
581  const int gg, const int bb) {
582  std::stringstream s;
583  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
584  THROW_MESSAGE(s.str());
586 }
587 
588 template <int Tensor_Dim0, int Tensor_Dim1>
590  Tensor_Dim0, Tensor_Dim1>
592  FieldApproximationBase base) {
593  std::stringstream s;
594  s << "Template for tensor dimension " << Tensor_Dim0 << "x" << Tensor_Dim1
595  << " not implemented";
596  THROW_MESSAGE(s.str());
598 }
599 
600 template <int Tensor_Dim0, int Tensor_Dim1>
602  Tensor_Dim0, Tensor_Dim1>
604  const int gg,
605  const int bb) {
606  std::stringstream s;
607  s << "Template for tensor dimension " << Tensor_Dim0 << "x" << Tensor_Dim1
608  << " not implemented";
609  THROW_MESSAGE(s.str());
611 }
612 
613 template <>
615 DataForcesAndSourcesCore::EntData::getFTensor1N<3>(
616  FieldApproximationBase base) {
617  double *t_n_ptr = &*getN(base).data().begin();
618  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(t_n_ptr, // HVEC0
619  &t_n_ptr[HVEC1],
620  &t_n_ptr[HVEC2]);
621 }
622 
623 template <>
625 DataForcesAndSourcesCore::EntData::getFTensor1N<3>(FieldApproximationBase base,
626  const int gg, const int bb) {
627  double *t_n_ptr = &getN(base)(gg, 3 * bb);
628  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(t_n_ptr, // HVEC0
629  &t_n_ptr[HVEC1],
630  &t_n_ptr[HVEC2]);
631 }
632 
633 template <>
635 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 3>(
636  FieldApproximationBase base) {
637  double *t_diff_n_ptr = &*getDiffN(base).data().begin();
639  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC0_2],
640  &t_diff_n_ptr[HVEC1_0], &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC1_2],
641  &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1], &t_diff_n_ptr[HVEC2_2]);
642 }
643 
644 template <>
646 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 3>(
647  FieldApproximationBase base, const int gg, const int bb) {
648  double *t_diff_n_ptr = &getDiffN(base)(gg, 9 * bb);
650  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC0_2],
651  &t_diff_n_ptr[HVEC1_0], &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC1_2],
652  &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1], &t_diff_n_ptr[HVEC2_2]);
653 }
654 
655 template <>
657 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 2>(
658  FieldApproximationBase base) {
659  double *t_diff_n_ptr = &*getDiffN(base).data().begin();
661  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC1_0],
662  &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1]);
663 }
664 
665 template <>
667 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 2>(
668  FieldApproximationBase base, const int gg, const int bb) {
669  double *t_diff_n_ptr = &getDiffN(base)(gg, 6 * bb);
671  t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC1_0],
672  &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1]);
673 }
674 
675 template <int Tensor_Dim0, int Tensor_Dim1>
677  Tensor_Dim0, Tensor_Dim1>
679  std::stringstream s;
680  s << "Template for tensor dimension " << Tensor_Dim0 << ", " << Tensor_Dim1
681  << " not implemented";
682  THROW_MESSAGE(s.str());
684  Tensor_Dim0, Tensor_Dim1>();
685 }
686 
687 template <>
689 DataForcesAndSourcesCore::EntData::getFTensor2N<3, 3>(
690  FieldApproximationBase base) {
691  double *t_n_ptr = &*(getN(base).data().begin());
693 
694  &t_n_ptr[HVEC0], &t_n_ptr[HVEC1], &t_n_ptr[HVEC2],
695 
696  &t_n_ptr[3 + HVEC0], &t_n_ptr[3 + HVEC1], &t_n_ptr[3 + HVEC2],
697 
698  &t_n_ptr[6 + HVEC0], &t_n_ptr[6 + HVEC1], &t_n_ptr[6 + HVEC2]
699 
700  );
701 }
702 
703 template <int Tensor_Dim0, int Tensor_Dim1>
705  Tensor_Dim0, Tensor_Dim1>
707  const int gg, const int bb) {
708  std::stringstream s;
709  s << "Template for tensor dimension " << Tensor_Dim0 << ", " << Tensor_Dim1
710  << " not implemented";
711  THROW_MESSAGE(s.str());
713  Tensor_Dim0, Tensor_Dim1>();
714 }
715 
716 template <>
718 DataForcesAndSourcesCore::EntData::getFTensor2N<3, 3>(
719  FieldApproximationBase base, const int gg, const int bb) {
720  double *t_n_ptr = &getN(base)(gg, 9 * bb);
722 
723  &t_n_ptr[HVEC0], &t_n_ptr[HVEC1], &t_n_ptr[HVEC2],
724 
725  &t_n_ptr[3 + HVEC0], &t_n_ptr[3 + HVEC1], &t_n_ptr[3 + HVEC2],
726 
727  &t_n_ptr[6 + HVEC0], &t_n_ptr[6 + HVEC1], &t_n_ptr[6 + HVEC2]
728 
729  );
730 }
731 
732 /**@}*/
733 
734 /** \name Bernstein-Bezier base only functions */
735 
736 /**@{*/
737 
738 boost::shared_ptr<MatrixInt> &
740  const std::string &field_name) {
741  return bbAlphaInduces[field_name];
742 }
743 
744 boost::shared_ptr<MatrixDouble> &
746  const std::string &field_name) {
747  return bbN[field_name];
748 }
749 
750 /**
751  * Get shared pointer to BB base base functions
752  */
753 const boost::shared_ptr<MatrixDouble> &
755  const std::string &field_name) const {
756  return bbN.at(field_name);
757 }
758 
759 /**
760  * Get shared pointer to BB derivatives of base base functions
761  */
762 boost::shared_ptr<MatrixDouble> &
764  const std::string &field_name) {
765  return bbDiffN[field_name];
766 }
767 
768 /**
769  * Get shared pointer to derivatives of BB base base functions
770  */
771 const boost::shared_ptr<MatrixDouble> &
773  const std::string &field_name) const {
774  return bbDiffN.at(field_name);
775 }
776 
777 std::map<std::string, boost::shared_ptr<MatrixInt>> &
779  return bbAlphaInduces;
780 }
781 
782 std::map<std::string, boost::shared_ptr<MatrixDouble>> &
784  return bbN;
785 }
786 
787 std::map<std::string, boost::shared_ptr<MatrixDouble>> &
789  return bbDiffN;
790 }
791 
792 boost::shared_ptr<MatrixInt> &
794  const size_t o) {
795  return bbAlphaInducesByOrder[o];
796 }
797 
798 boost::shared_ptr<MatrixDouble> &
800  return bbNByOrder[o];
801 }
802 
803 const boost::shared_ptr<MatrixDouble> &
805  return bbNByOrder[o];
806 }
807 
808 boost::shared_ptr<MatrixDouble> &
810  return bbDiffNByOrder[o];
811 }
812 
813 const boost::shared_ptr<MatrixDouble> &
815  const size_t o) const {
816  return bbDiffNByOrder[o];
817 }
818 
819 std::array<boost::shared_ptr<MatrixInt>,
822  return bbAlphaInducesByOrder;
823 }
824 
825 std::array<boost::shared_ptr<MatrixDouble>,
828  return bbNByOrder;
829 }
830 
831 std::array<boost::shared_ptr<MatrixDouble>,
834  return bbDiffNByOrder;
835 }
836 
837 boost::shared_ptr<MatrixInt> &
839  const std::string &field_name) {
840  return entDataPtr->getBBAlphaIndicesSharedPtr(field_name);
841 }
842 
843 boost::shared_ptr<MatrixDouble> &
845  const std::string &field_name) {
846  return entDataPtr->getBBNSharedPtr(field_name);
847 }
848 
849 const boost::shared_ptr<MatrixDouble> &
851  const std::string &field_name) const {
852  return entDataPtr->getBBNSharedPtr(field_name);
853 }
854 
855 /**
856  * Get shared pointer to BB derivatives of base base functions
857  */
858 boost::shared_ptr<MatrixDouble> &
860  const std::string &field_name) {
861  return entDataPtr->getBBDiffNSharedPtr(field_name);
862 }
863 
864 /**
865  * Get shared pointer to derivatives of BB base base functions
866  */
867 const boost::shared_ptr<MatrixDouble> &
869  const std::string &field_name) const {
870  return entDataPtr->getBBDiffNSharedPtr(field_name);
871 }
872 
873 /**@}*/
874 
875 
876 } // namespace MoFEM
877 
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:501
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:477
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:620
MoFEMErrorCode setElementType(const EntityType type)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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:144
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:596
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:407
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.