11 : sEnse(0), oRder(0), bAse(
NOBASE), entDataBitRefLevel(),
12 N(baseFunctionsAndBaseDerivatives[ZeroDerivative]),
13 diffN(baseFunctionsAndBaseDerivatives[FirstDerivative]) {
14 if (allocate_base_matrices) {
17 for (
int b = 0; b !=
LASTBASE; ++b) {
26boost::shared_ptr<MatrixDouble> &
29 return baseFunctionsAndBaseDerivatives[direvatie][base];
32boost::shared_ptr<MatrixDouble> &
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);
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;
58 for (
auto tt = moab::CN::TypeDimensionMap[dd].first;
59 tt <= moab::CN::TypeDimensionMap[dd].second; ++tt) {
87 const boost::shared_ptr<EntitiesFieldData> &data_ptr) {
92 for (
int tt = MBVERTEX; tt != MBMAXTYPE; ++tt) {
93 auto &ent_data = data_ptr->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));
99 derived_ent_data.resize(ent_data.size());
104 const boost::shared_ptr<EntitiesFieldData> &data_ptr)
122 dOfs.resize(0,
false);
131 CHKERR e.resetFieldDependentData();
139 auto make_swap = [](boost::shared_ptr<MatrixDouble> &ptr,
140 boost::shared_ptr<MatrixDouble> &ptrBB,
141 boost::shared_ptr<MatrixDouble> &swap_ptr) {
150 make_swap(getNSharedPtr(base), getBBNSharedPtr(
field_name), swapBaseNPtr);
151 make_swap(getDiffNSharedPtr(base), getBBDiffNSharedPtr(
field_name),
160 for (
int tt = MBVERTEX; tt != MBENTITYSET; ++tt) {
162 for (
auto &side_data : ent_data)
171 auto make_swap = [](boost::shared_ptr<MatrixDouble> &ptr,
172 boost::shared_ptr<MatrixDouble> &ptrBB,
173 boost::shared_ptr<MatrixDouble> &swap_ptr) {
182 make_swap(getDerivedNSharedPtr(base), getBBNSharedPtr(
field_name),
184 make_swap(getDerivedDiffNSharedPtr(base), getBBDiffNSharedPtr(
field_name),
190 const boost::shared_ptr<EntitiesFieldData::EntData> &ent_data_ptr)
194 return entDataPtr->getSense();
197boost::shared_ptr<MatrixDouble> &
200 if (baseFunctionsAndBaseDerivatives[derivative][base])
201 return baseFunctionsAndBaseDerivatives[derivative][base];
203 return entDataPtr->getNSharedPtr(base, derivative);
206boost::shared_ptr<MatrixDouble> &
212 return entDataPtr->getNSharedPtr(base);
215boost::shared_ptr<MatrixDouble> &
221 return entDataPtr->getDiffNSharedPtr(base);
223const boost::shared_ptr<MatrixDouble> &
229 return entDataPtr->getNSharedPtr(base);
231const boost::shared_ptr<MatrixDouble> &
237 return entDataPtr->getDiffNSharedPtr(base);
242 os <<
"sEnse: " << e.
getSense() << std::endl
243 <<
"oRder: " << e.
getOrder() << std::endl
244 <<
"global indices: " << e.
getIndices() << std::endl
247 os <<
"fieldData: " << std::fixed << std::setprecision(2) << e.
getFieldData()
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)
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;
265 os <<
"N: " << std::fixed << base << std::endl
266 <<
"diffN: " << std::fixed << diff_base;
273 os <<
"dataOnEntities[" << moab::CN::EntityTypeName(
t) <<
"][" << nn
287EntitiesFieldData::EntData::getFTensor1FieldData<3>() {
288 if (dOfs[0]->getNbOfCoeffs() != 3) {
290 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
291 s <<
" but you ask for tensor rank 1 dimension 3";
294 double *ptr = &*fieldData.data().begin();
301EntitiesFieldData::EntData::getFTensor1FieldData<2>() {
302 if (dOfs[0]->getNbOfCoeffs() != 2) {
304 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
305 s <<
" but you ask for tensor rank 1 dimension 2";
308 double *ptr = &*fieldData.data().begin();
314EntitiesFieldData::EntData::getFTensor1FieldData<1>() {
315 if (dOfs[0]->getNbOfCoeffs() != 1) {
317 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
318 s <<
" but you ask for tensor rank 1 dimension 1";
321 double *ptr = &*fieldData.data().begin();
327EntitiesFieldData::EntData::getFTensor2FieldData<1, 1>() {
328 if (dOfs[0]->getNbOfCoeffs() != 1) {
330 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
331 s <<
" but you ask for tensor rank 2 dimensions 1 by 1 so 1 coefficients "
335 double *ptr = &*fieldData.data().begin();
341EntitiesFieldData::EntData::getFTensor2FieldData<1, 2>() {
342 if (dOfs[0]->getNbOfCoeffs() != 1) {
344 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
345 s <<
" but you ask for tensor rank 2 dimensions 1 by 2 so 2 coefficients "
349 double *ptr = &*fieldData.data().begin();
355EntitiesFieldData::EntData::getFTensor2FieldData<1, 3>() {
356 if (dOfs[0]->getNbOfCoeffs() != 1) {
358 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
359 s <<
" but you ask for tensor rank 2 dimensions 1 by 3 so 3 coefficients "
363 double *ptr = &*fieldData.data().begin();
370EntitiesFieldData::EntData::getFTensor2FieldData<2, 2>() {
371 if (dOfs[0]->getNbOfCoeffs() != 4) {
373 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
374 s <<
" but you ask for tensor rank 2 dimensions 2 by 2 so 4 coefficients "
378 double *ptr = &*fieldData.data().begin();
380 ptr, &ptr[1], &ptr[2], &ptr[3]);
385EntitiesFieldData::EntData::getFTensor2FieldData<3, 3>() {
386 if (dOfs[0]->getNbOfCoeffs() != 9) {
388 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
389 s <<
" but you ask for tensor rank 2 dimensions 3 by 3 so 9 coefficients "
393 double *ptr = &*fieldData.data().begin();
395 ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
401EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<3>() {
402 if (dOfs[0]->getNbOfCoeffs() != 6) {
404 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
405 s <<
" but you ask for symmetric tensor rank 2 dimensions 3 by 3 so 6 "
410 double *ptr = &*fieldData.data().begin();
412 ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
417EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<2>() {
418 if (dOfs[0]->getNbOfCoeffs() != 3) {
420 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
421 s <<
" but you ask for symmetric tensor rank 2 dimensions 2 by 2 so 3 "
426 double *ptr = &*fieldData.data().begin();
428 ptr, &ptr[1], &ptr[2]);
433 if (dOfs[0]->getNbOfCoeffs() != 1) {
435 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
436 s <<
" but expected scalar field, tensor of rank 0";
440 &*fieldData.data().begin());
443template <
int Tensor_Dim>
448 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
453template <
int Tensor_Dim>
456 return getFTensor1DiffN<Tensor_Dim>(
bAse);
464EntitiesFieldData::EntData::getFTensor1DiffN<3>(
466 double *ptr = &*getDiffN(base).data().begin();
476EntitiesFieldData::EntData::getFTensor1DiffN<3>() {
477 return getFTensor1DiffN<3>(
bAse);
485EntitiesFieldData::EntData::getFTensor1DiffN<2>(
487 double *ptr = &*getDiffN(base).data().begin();
496EntitiesFieldData::EntData::getFTensor1DiffN<2>() {
497 return getFTensor1DiffN<2>(
bAse);
500template <
int Tensor_Dim>
503 const int gg,
const int bb) {
505 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
515EntitiesFieldData::EntData::getFTensor1DiffN<3>(
517 double *ptr = &getDiffN(base)(gg, 3 * bb);
527EntitiesFieldData::EntData::getFTensor1DiffN<3>(
const int gg,
const int bb) {
528 return getFTensor1DiffN<3>(
bAse, gg, bb);
536EntitiesFieldData::EntData::getFTensor1DiffN<2>(
538 double *ptr = &getDiffN(base)(gg, 2 * bb);
547EntitiesFieldData::EntData::getFTensor1DiffN<2>(
const int gg,
const int bb) {
548 return getFTensor1DiffN<2>(
bAse, gg, bb);
557template <
int Tensor_Dim>
561 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
566template <
int Tensor_Dim>
569 const int gg,
const int bb) {
571 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
576template <
int Tensor_Dim0,
int Tensor_Dim1>
578 Tensor_Dim0, Tensor_Dim1>
581 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
"x" << Tensor_Dim1
582 <<
" not implemented";
587template <
int Tensor_Dim0,
int Tensor_Dim1>
589 Tensor_Dim0, Tensor_Dim1>
591 const int gg,
const int bb) {
593 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
"x" << Tensor_Dim1
594 <<
" not implemented";
602 double *t_n_ptr = &*getN(base).
data().begin();
611 const int gg,
const int bb) {
612 double *t_n_ptr = &getN(base)(gg, 3 * bb);
620EntitiesFieldData::EntData::getFTensor2DiffN<3, 3>(
622 double *t_diff_n_ptr = &*getDiffN(base).
data().begin();
632 const int gg,
const int bb) {
633 double *t_diff_n_ptr = &getDiffN(base)(gg, 9 * bb);
642EntitiesFieldData::EntData::getFTensor2DiffN<3, 2>(
644 double *t_diff_n_ptr = &*getDiffN(base).data().begin();
653 const int gg,
const int bb) {
654 double *t_diff_n_ptr = &getDiffN(base)(gg, 6 * bb);
663 double *ptr = &(getN(base, BaseDerivatives::SecondDerivative))(0, 0);
678template <
int Tensor_Dim0,
int Tensor_Dim1>
680 Tensor_Dim0, Tensor_Dim1>
683 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
", " << Tensor_Dim1
684 <<
" not implemented";
687 Tensor_Dim0, Tensor_Dim1>();
693 double *t_n_ptr = &*(getN(base).data().begin());
705template <
int Tensor_Dim0,
int Tensor_Dim1>
707 Tensor_Dim0, Tensor_Dim1>
709 const int gg,
const int bb) {
711 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
", " << Tensor_Dim1
712 <<
" not implemented";
715 Tensor_Dim0, Tensor_Dim1>();
721 const int gg,
const int bb) {
722 double *t_n_ptr = &getN(base)(gg, 9 * bb);
740boost::shared_ptr<MatrixInt> &
746boost::shared_ptr<MatrixDouble> &
754const boost::shared_ptr<MatrixDouble> &
763boost::shared_ptr<MatrixDouble> &
771const boost::shared_ptr<MatrixDouble> &
777std::map<std::string, boost::shared_ptr<MatrixInt>> &
779 return bbAlphaIndices;
782std::map<std::string, boost::shared_ptr<MatrixDouble>> &
787std::map<std::string, boost::shared_ptr<MatrixDouble>> &
792boost::shared_ptr<MatrixInt> &
794 return bbAlphaIndicesByOrder[
o];
797boost::shared_ptr<MatrixDouble> &
799 return bbNByOrder[
o];
802boost::shared_ptr<MatrixDouble> &
804 return bbDiffNByOrder[
o];
807std::array<boost::shared_ptr<MatrixInt>,
810 return bbAlphaIndicesByOrder;
813std::array<boost::shared_ptr<MatrixDouble>,
819std::array<boost::shared_ptr<MatrixDouble>,
822 return bbDiffNByOrder;
825boost::shared_ptr<MatrixInt> &
828 return entDataPtr->getBBAlphaIndicesSharedPtr(
field_name);
831boost::shared_ptr<MatrixDouble> &
834 return entDataPtr->getBBNSharedPtr(
field_name);
837const boost::shared_ptr<MatrixDouble> &
840 return entDataPtr->getBBNSharedPtr(
field_name);
846boost::shared_ptr<MatrixDouble> &
849 return entDataPtr->getBBDiffNSharedPtr(
field_name);
855const boost::shared_ptr<MatrixDouble> &
858 return entDataPtr->getBBDiffNSharedPtr(
field_name);
864 return entDataBitRefLevel;
867std::vector<BitRefLevel> &
869 return entDataPtr->getEntDataBitRefLevel();
T data[Tensor_Dim0][Tensor_Dim1]
FieldApproximationBase
approximation base
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
const double c
speed of light (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
static void constructor_derived_data(DerivedEntitiesFieldData *derived_data, const boost::shared_ptr< EntitiesFieldData > &data_ptr)
static void constructor_data(EntitiesFieldData *data, const EntityType type)
constexpr double t
plate stiffness
constexpr auto field_name
Derived ata on single entity (This is passed as argument to DataOperator::doWork)
boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Used by Bernstein base to keep temporally pointer.
std::vector< BitRefLevel > & getEntDataBitRefLevel()
int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
DerivedEntData(const boost::shared_ptr< EntitiesFieldData::EntData > &ent_data_ptr)
boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives derivative)
this class derive data form other data structure
MoFEMErrorCode setElementType(const EntityType type)
const boost::shared_ptr< EntitiesFieldData > dataPtr
DerivedEntitiesFieldData(const boost::shared_ptr< EntitiesFieldData > &data_ptr)
Data on single entity (This is passed as argument to DataOperator::doWork)
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
FieldSpace sPace
Entity space.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBNByOrderArray()
auto getFTensor2N()
Get base functions for Hdiv space.
VectorDouble fieldData
Field data on entity.
ApproximationOrder getOrder() const
get approximation order
VectorInt localIndices
Local indices on entity.
virtual std::array< boost::shared_ptr< MatrixInt >, MaxBernsteinBezierOrder > & getBBAlphaIndicesByOrderArray()
virtual std::map< std::string, boost::shared_ptr< MatrixInt > > & getBBAlphaIndicesMap()
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN()
Get derivatives of base functions.
static constexpr size_t MaxBernsteinBezierOrder
VectorInt iNdices
Global indices on entity.
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesSharedPtr(const std::string &field_name)
const VectorDouble & getFieldData() const
get dofs values
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNByOrderSharedPtr(const size_t o)
get BB base derivative by order
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBNMap()
get hash map of base function for BB base, key is a field name
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Swap bases functions.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
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
virtual std::vector< BitRefLevel > & getEntDataBitRefLevel()
virtual int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
VectorFieldEntities fieldEntities
Field entities.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
virtual boost::shared_ptr< MatrixInt > & getBBAlphaIndicesByOrderSharedPtr(const size_t o)
get ALpha indices for BB base by order
MoFEMErrorCode resetFieldDependentData()
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.
auto getFTensor1N()
Get base functions for Hdiv space.
virtual boost::shared_ptr< MatrixDouble > & getBBNSharedPtr(const std::string &field_name)
std::array< std::array< boost::shared_ptr< MatrixDouble >, LASTBASE >, LastDerivative > baseFunctionsAndBaseDerivatives
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN()
Get derivatives of base functions for Hdiv space.
VectorDofs dOfs
DoFs on entity.
virtual boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
virtual boost::shared_ptr< MatrixDouble > & getBBNByOrderSharedPtr(const size_t o)
get BB base by order
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives direvatie)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
EntData(const bool allocate_base_matrices=true)
FieldApproximationBase bAse
Field approximation base.
virtual boost::shared_ptr< MatrixDouble > & getBBDiffNSharedPtr(const std::string &field_name)
data structure for finite element entity
friend std::ostream & operator<<(std::ostream &os, const EntitiesFieldData &e)
MoFEMErrorCode resetFieldDependentData()
std::bitset< LASTBASE > bAse
bases on element
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
virtual MoFEMErrorCode baseSwap(const std::string &field_name, const FieldApproximationBase base)
Swap approximation base.
virtual MoFEMErrorCode setElementType(const EntityType type)