23 : sEnse(0), oRder(0), bAse(
NOBASE) {
24 if (allocate_base_matrices)
33 boost::shared_ptr<MatrixDouble> &
39 const boost::shared_ptr<MatrixDouble> &
45 boost::shared_ptr<MatrixDouble> &
51 const boost::shared_ptr<MatrixDouble> &
58 const EntityType
type) {
69 for (
int ee = 0; ee != 6; ++ee) {
72 for (
int ff = 0; ff != 4; ++ff) {
79 for (
int ee = 0; ee != 3; ++ee) {
86 for (
int ee = 0; ee != 4; ++ee) {
100 for (
int ee = 0; ee != 9; ++ee) {
103 for (
int qq = 0; qq != 5; ++qq) {
106 for (
int ff = 0; ff != 5; ++ff) {
130 const boost::shared_ptr<DataForcesAndSourcesCore> &data_ptr) {
135 for (
int tt = MBVERTEX; tt != MBMAXTYPE; ++tt) {
136 auto &ent_data = data_ptr->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));
146 const boost::shared_ptr<DataForcesAndSourcesCore> &data_ptr)
154 for (EntityType tt = MBVERTEX; tt != MBMAXTYPE; ++tt)
166 dOfs.resize(0,
false);
173 for (EntityType t = MBVERTEX; t != MBMAXTYPE; t++)
175 CHKERR e.resetFieldDependentData();
183 auto make_swap = [](boost::shared_ptr<MatrixDouble> &ptr,
184 boost::shared_ptr<MatrixDouble> &ptrBB,
185 boost::shared_ptr<MatrixDouble> &swap_ptr) {
194 make_swap(getNSharedPtr(base), getBBNSharedPtr(field_name), swapBaseNPtr);
195 make_swap(getDiffNSharedPtr(base), getBBDiffNSharedPtr(field_name),
203 for (
int tt = MBVERTEX; tt != MBMAXTYPE; ++tt) {
205 for (
auto &side_data : ent_data)
206 CHKERR side_data.baseSwap(field_name, base);
214 auto make_swap = [](boost::shared_ptr<MatrixDouble> &ptr,
215 boost::shared_ptr<MatrixDouble> &ptrBB,
216 boost::shared_ptr<MatrixDouble> &swap_ptr) {
225 make_swap(getDerivedNSharedPtr(base), getBBNSharedPtr(field_name),
227 make_swap(getDerivedDiffNSharedPtr(base), getBBDiffNSharedPtr(field_name),
233 const boost::shared_ptr<DataForcesAndSourcesCore::EntData> &ent_data_ptr)
237 return entDataPtr->getSense();
240 boost::shared_ptr<MatrixDouble> &
246 return entDataPtr->getNSharedPtr(base);
248 boost::shared_ptr<MatrixDouble> &
254 return entDataPtr->getDiffNSharedPtr(base);
256 const boost::shared_ptr<MatrixDouble> &
262 return entDataPtr->getNSharedPtr(base);
264 const boost::shared_ptr<MatrixDouble> &
270 return entDataPtr->getDiffNSharedPtr(base);
275 os <<
"sEnse: " << e.
getSense() << std::endl
276 <<
"oRder: " << e.
getOrder() << std::endl
277 <<
"global indices: " << e.
getIndices() << std::endl
280 os <<
"fieldData: " << std::fixed << std::setprecision(2) << e.
getFieldData() << std::endl;
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)
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;
297 os <<
"N: " << std::fixed << base << std::endl
298 <<
"diffN: " << std::fixed << diff_base;
303 for (
unsigned int nn = 0; nn < e.
dataOnEntities[MBVERTEX].size(); nn++) {
304 os <<
"dataOnEntities[MBVERTEX][" << nn <<
"]" << std::endl
307 for (
unsigned int ee = 0; ee < e.
dataOnEntities[MBEDGE].size(); ee++) {
308 os <<
"dataOnEntities[MBEDGE][" << ee <<
"]" << std::endl
311 for (
unsigned int ff = 0; ff < e.
dataOnEntities[MBTRI].size(); ff++) {
312 os <<
"dataOnEntities[MBTRI][" << ff <<
"] " << std::endl
315 for (
unsigned int vv = 0; vv < e.
dataOnEntities[MBTET].size(); vv++) {
316 os <<
"dataOnEntities[MBTET][" << vv <<
"]" << std::endl
328 DataForcesAndSourcesCore::EntData::getFTensor1FieldData<3>() {
329 if (dOfs[0]->getNbOfCoeffs() != 3) {
331 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
332 s <<
" but you ask for tensor rank 1 dimension 3";
335 double *ptr = &*fieldData.data().begin();
342 DataForcesAndSourcesCore::EntData::getFTensor1FieldData<2>() {
343 if (dOfs[0]->getNbOfCoeffs() != 2) {
345 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
346 s <<
" but you ask for tensor rank 1 dimension 3";
349 double *ptr = &*fieldData.data().begin();
355 DataForcesAndSourcesCore::EntData::getFTensor2FieldData<3, 3>() {
356 if (dOfs[0]->getNbOfCoeffs() != 9) {
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 "
363 double *ptr = &*fieldData.data().begin();
365 ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
371 DataForcesAndSourcesCore::EntData::getFTensor2SymmetricFieldData<3>() {
372 if (dOfs[0]->getNbOfCoeffs() != 6) {
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 "
380 double *ptr = &*fieldData.data().begin();
382 ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
387 DataForcesAndSourcesCore::EntData::getFTensor2SymmetricFieldData<2>() {
388 if (dOfs[0]->getNbOfCoeffs() != 3) {
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 "
396 double *ptr = &*fieldData.data().begin();
398 ptr, &ptr[1], &ptr[2]);
403 if (dOfs[0]->getNbOfCoeffs() != 1) {
405 s <<
"Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
406 s <<
" but expected scalar field, tensor of rank 0";
410 &*fieldData.data().begin());
413 template <
int Tensor_Dim>
418 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
423 template <
int Tensor_Dim>
426 return getFTensor1DiffN<Tensor_Dim>(
bAse);
429 template <
int Tensor_Dim>
434 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
439 template <
int Tensor_Dim>
442 return getFTensor1DiffN<Tensor_Dim>(
bAse, bb);
450 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
452 double *ptr = &*getDiffN(base).
data().begin();
461 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>() {
462 return getFTensor1DiffN<3>(
bAse);
470 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
472 double *ptr = &*getDiffN(base).
data().begin();
474 &ptr[3 * bb], &ptr[3 * bb + 1], &ptr[3 * bb + 2], getDiffN(base).size2());
482 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
const int bb) {
483 return getFTensor1DiffN<3>(
bAse, bb);
491 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
493 double *ptr = &*getDiffN(base).
data().begin();
502 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>() {
503 return getFTensor1DiffN<2>(
bAse);
511 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
513 double *ptr = &*getDiffN(base).
data().begin();
515 getDiffN(base).size1());
523 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
const int bb) {
524 return getFTensor1DiffN<2>(
bAse, bb);
527 template <
int Tensor_Dim>
532 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
542 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
544 double *ptr = &getDiffN(base)(gg, 3 * bb);
553 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
const int gg,
555 return getFTensor1DiffN<3>(
bAse, gg, bb);
563 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
565 double *ptr = &getDiffN(base)(gg, 2 * bb);
574 DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
const int gg,
576 return getFTensor1DiffN<2>(
bAse, gg, bb);
585 template <
int Tensor_Dim>
589 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
594 template <
int Tensor_Dim>
597 const int gg,
const int bb) {
599 s <<
"Template for tensor dimension " << Tensor_Dim <<
" not implemented";
604 template <
int Tensor_Dim0,
int Tensor_Dim1>
606 Tensor_Dim0, Tensor_Dim1>
610 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
"x" << Tensor_Dim1
611 <<
" not implemented";
616 template <
int Tensor_Dim0,
int Tensor_Dim1>
618 Tensor_Dim0, Tensor_Dim1>
623 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
"x" << Tensor_Dim1
624 <<
" not implemented";
631 DataForcesAndSourcesCore::EntData::getFTensor1N<3>(
633 double *t_n_ptr = &*getN(base).
data().begin();
642 const int gg,
const int bb) {
643 double *t_n_ptr = &getN(base)(gg, 3 * bb);
651 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 3>(
653 double *t_diff_n_ptr = &*getDiffN(base).
data().begin();
662 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 3>(
664 double *t_diff_n_ptr = &getDiffN(base)(gg, 9 * bb);
673 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 2>(
675 double *t_diff_n_ptr = &*getDiffN(base).
data().begin();
683 DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 2>(
685 double *t_diff_n_ptr = &getDiffN(base)(gg, 6 * bb);
691 template <
int Tensor_Dim0,
int Tensor_Dim1>
693 Tensor_Dim0, Tensor_Dim1>
696 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
", " << Tensor_Dim1
697 <<
" not implemented";
700 Tensor_Dim0, Tensor_Dim1>();
705 DataForcesAndSourcesCore::EntData::getFTensor2N<3, 3>(
707 double *t_n_ptr = &*(getN(base).data().begin());
719 template <
int Tensor_Dim0,
int Tensor_Dim1>
721 Tensor_Dim0, Tensor_Dim1>
723 const int gg,
const int bb) {
725 s <<
"Template for tensor dimension " << Tensor_Dim0 <<
", " << Tensor_Dim1
726 <<
" not implemented";
729 Tensor_Dim0, Tensor_Dim1>();
734 DataForcesAndSourcesCore::EntData::getFTensor2N<3, 3>(
736 double *t_n_ptr = &getN(base)(gg, 9 * bb);
754 boost::shared_ptr<MatrixInt> &
756 const std::string &field_name) {
757 return bbAlphaInduces[field_name];
760 boost::shared_ptr<MatrixDouble> &
762 const std::string &field_name) {
763 return bbN[field_name];
769 const boost::shared_ptr<MatrixDouble> &
771 const std::string &field_name)
const {
772 return bbN.at(field_name);
778 boost::shared_ptr<MatrixDouble> &
780 const std::string &field_name) {
781 return bbDiffN[field_name];
787 const boost::shared_ptr<MatrixDouble> &
789 const std::string &field_name)
const {
790 return bbDiffN.at(field_name);
793 std::map<std::string, boost::shared_ptr<MatrixInt>> &
795 return bbAlphaInduces;
798 std::map<std::string, boost::shared_ptr<MatrixDouble>> &
803 std::map<std::string, boost::shared_ptr<MatrixDouble>> &
808 boost::shared_ptr<MatrixInt> &
811 return bbAlphaInducesByOrder[
o];
814 boost::shared_ptr<MatrixDouble> &
816 return bbNByOrder[
o];
819 const boost::shared_ptr<MatrixDouble> &
821 return bbNByOrder[
o];
824 boost::shared_ptr<MatrixDouble> &
826 return bbDiffNByOrder[
o];
829 const boost::shared_ptr<MatrixDouble> &
831 const size_t o)
const {
832 return bbDiffNByOrder[
o];
835 std::array<boost::shared_ptr<MatrixInt>,
838 return bbAlphaInducesByOrder;
841 std::array<boost::shared_ptr<MatrixDouble>,
847 std::array<boost::shared_ptr<MatrixDouble>,
850 return bbDiffNByOrder;
853 boost::shared_ptr<MatrixInt> &
855 const std::string &field_name) {
856 return entDataPtr->getBBAlphaIndicesSharedPtr(field_name);
859 boost::shared_ptr<MatrixDouble> &
861 const std::string &field_name) {
862 return entDataPtr->getBBNSharedPtr(field_name);
865 const boost::shared_ptr<MatrixDouble> &
867 const std::string &field_name)
const {
868 return entDataPtr->getBBNSharedPtr(field_name);
874 boost::shared_ptr<MatrixDouble> &
876 const std::string &field_name) {
877 return entDataPtr->getBBDiffNSharedPtr(field_name);
883 const boost::shared_ptr<MatrixDouble> &
885 const std::string &field_name)
const {
886 return entDataPtr->getBBDiffNSharedPtr(field_name);
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(a)
Throw MoFEM exception.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
implementation of Data Operators for Forces and Sources
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
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
MoFEMErrorCode resetFieldDependentData()
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)
FieldSpace sPace
Entity space.
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
VectorDofs dOfs
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
DataForcesAndSourcesCore()
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)