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<3, 3>() {
328 if (dOfs[0]->getNbOfCoeffs() != 9) {
330 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
331 s <<
" but you ask for tensor rank 2 dimensions 3 by 3 so 9 coefficients "
335 double *ptr = &*fieldData.data().begin();
337 ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
343EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<3>() {
344 if (dOfs[0]->getNbOfCoeffs() != 6) {
346 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
347 s <<
" but you ask for symmetric tensor rank 2 dimensions 3 by 3 so 6 "
352 double *ptr = &*fieldData.data().begin();
354 ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
359EntitiesFieldData::EntData::getFTensor2SymmetricFieldData<2>() {
360 if (dOfs[0]->getNbOfCoeffs() != 3) {
362 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
363 s <<
" but you ask for symmetric tensor rank 2 dimensions 2 by 2 so 3 "
368 double *ptr = &*fieldData.data().begin();
370 ptr, &ptr[1], &ptr[2]);
375 if (dOfs[0]->getNbOfCoeffs() != 1) {
377 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
378 s <<
" but expected scalar field, tensor of rank 0";
382 &*fieldData.data().begin());
385template <
int Tensor_Dim>
390 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
395template <
int Tensor_Dim>
398 return getFTensor1DiffN<Tensor_Dim>(
bAse);
406EntitiesFieldData::EntData::getFTensor1DiffN<3>(
408 double *ptr = &*getDiffN(base).data().begin();
418EntitiesFieldData::EntData::getFTensor1DiffN<3>() {
419 return getFTensor1DiffN<3>(
bAse);
427EntitiesFieldData::EntData::getFTensor1DiffN<2>(
429 double *ptr = &*getDiffN(base).data().begin();
438EntitiesFieldData::EntData::getFTensor1DiffN<2>() {
439 return getFTensor1DiffN<2>(
bAse);
442template <
int Tensor_Dim>
445 const int gg,
const int bb) {
447 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
457EntitiesFieldData::EntData::getFTensor1DiffN<3>(
459 double *ptr = &getDiffN(base)(gg, 3 * bb);
469EntitiesFieldData::EntData::getFTensor1DiffN<3>(
const int gg,
const int bb) {
470 return getFTensor1DiffN<3>(
bAse, gg, bb);
478EntitiesFieldData::EntData::getFTensor1DiffN<2>(
480 double *ptr = &getDiffN(base)(gg, 2 * bb);
489EntitiesFieldData::EntData::getFTensor1DiffN<2>(
const int gg,
const int bb) {
490 return getFTensor1DiffN<2>(
bAse, gg, bb);
499template <
int Tensor_Dim>
503 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
508template <
int Tensor_Dim>
511 const int gg,
const int bb) {
513 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
518template <
int Tensor_Dim0,
int Tensor_Dim1>
520 Tensor_Dim0, Tensor_Dim1>
523 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
"x" << Tensor_Dim1
524 <<
" not implemented";
529template <
int Tensor_Dim0,
int Tensor_Dim1>
531 Tensor_Dim0, Tensor_Dim1>
533 const int gg,
const int bb) {
535 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
"x" << Tensor_Dim1
536 <<
" not implemented";
544 double *t_n_ptr = &*getN(base).
data().begin();
553 const int gg,
const int bb) {
554 double *t_n_ptr = &getN(base)(gg, 3 * bb);
562EntitiesFieldData::EntData::getFTensor2DiffN<3, 3>(
564 double *t_diff_n_ptr = &*getDiffN(base).
data().begin();
574 const int gg,
const int bb) {
575 double *t_diff_n_ptr = &getDiffN(base)(gg, 9 * bb);
584EntitiesFieldData::EntData::getFTensor2DiffN<3, 2>(
586 double *t_diff_n_ptr = &*getDiffN(base).data().begin();
595 const int gg,
const int bb) {
596 double *t_diff_n_ptr = &getDiffN(base)(gg, 6 * bb);
605 double *ptr = &(getN(base, BaseDerivatives::SecondDerivative))(0, 0);
620template <
int Tensor_Dim0,
int Tensor_Dim1>
622 Tensor_Dim0, Tensor_Dim1>
625 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
", " << Tensor_Dim1
626 <<
" not implemented";
629 Tensor_Dim0, Tensor_Dim1>();
635 double *t_n_ptr = &*(getN(base).data().begin());
647template <
int Tensor_Dim0,
int Tensor_Dim1>
649 Tensor_Dim0, Tensor_Dim1>
651 const int gg,
const int bb) {
653 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
", " << Tensor_Dim1
654 <<
" not implemented";
657 Tensor_Dim0, Tensor_Dim1>();
663 const int gg,
const int bb) {
664 double *t_n_ptr = &getN(base)(gg, 9 * bb);
682boost::shared_ptr<MatrixInt> &
688boost::shared_ptr<MatrixDouble> &
696const boost::shared_ptr<MatrixDouble> &
705boost::shared_ptr<MatrixDouble> &
713const boost::shared_ptr<MatrixDouble> &
719std::map<std::string, boost::shared_ptr<MatrixInt>> &
721 return bbAlphaIndices;
724std::map<std::string, boost::shared_ptr<MatrixDouble>> &
729std::map<std::string, boost::shared_ptr<MatrixDouble>> &
734boost::shared_ptr<MatrixInt> &
736 return bbAlphaIndicesByOrder[
o];
739boost::shared_ptr<MatrixDouble> &
741 return bbNByOrder[
o];
744boost::shared_ptr<MatrixDouble> &
746 return bbDiffNByOrder[
o];
749std::array<boost::shared_ptr<MatrixInt>,
752 return bbAlphaIndicesByOrder;
755std::array<boost::shared_ptr<MatrixDouble>,
761std::array<boost::shared_ptr<MatrixDouble>,
764 return bbDiffNByOrder;
767boost::shared_ptr<MatrixInt> &
770 return entDataPtr->getBBAlphaIndicesSharedPtr(
field_name);
773boost::shared_ptr<MatrixDouble> &
776 return entDataPtr->getBBNSharedPtr(
field_name);
779const boost::shared_ptr<MatrixDouble> &
782 return entDataPtr->getBBNSharedPtr(
field_name);
788boost::shared_ptr<MatrixDouble> &
791 return entDataPtr->getBBDiffNSharedPtr(
field_name);
797const boost::shared_ptr<MatrixDouble> &
800 return entDataPtr->getBBDiffNSharedPtr(
field_name);
806 return entDataBitRefLevel;
809std::vector<BitRefLevel> &
811 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)