10 boost::shared_ptr<MatrixDouble> jac_ptr)
14 for (
auto t = MBEDGE;
t != MBMAXTYPE; ++
t)
26 const auto nb_base_functions = data.
getN(
NOBASE).size2();
27 if (nb_base_functions) {
34 if (nb_gauss_pts != data.
getDiffN().size1())
36 "Inconsistent number base functions and gauss points");
37 if (nb_base_functions != data.
getDiffN().size2() / 3)
39 "Inconsistent number of base functions");
40 if (coords.size() != 3 * nb_base_functions)
42 "Number of vertex coordinates and base functions is inconsistent "
44 coords.size() / 3, nb_base_functions);
47 jacPtr->resize(9, nb_gauss_pts,
false);
56 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
58 &coords[0], &coords[1], &coords[2]);
59 for (
size_t bb = 0; bb != nb_base_functions; ++bb) {
60 t_jac(
i,
j) += t_coords(
i) * (t_vol_inv_jac(
k,
j) * t_diff_base(
k));
90 CHKERR transform_base(*(
m.second));
95 CHKERR transform_base(*ptr);
119 data.
getDiffN(base).size2(),
false);
121 unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
122 unsigned int nb_base_functions = data.
getDiffN(base).size2() / 9;
123 if (nb_base_functions == 0)
129 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
135 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
136 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
137 t_inv_diff_n(
i,
j) = t_inv_jac(
k,
j) * t_diff_n(
i,
k);
167 for (
size_t gg = 0; gg != nb_int_pts; ++gg) {
168 t_w *= sqrt(t_normal(
i) * t_normal(
i)) /
a;
174 "Number of rows in getNormalsAtGaussPts should be equal to "
175 "number of integration points, but is not, i.e. %ld != %ld",
192 for (
size_t gg = 0; gg != nb_int_pts; ++gg) {
193 t_w *= sqrt(t_tangent(
i) * t_tangent(
i)) /
a;
199 "Number of rows in getTangentAtGaussPts should be equal to "
200 "number of integration points, but is not, i.e. %ld != %ld",
211 const auto nb_integration_pts =
detPtr->size();
216 "Inconsistent number of data points");
221 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
242 "Pointer for detPtr not allocated");
245 "Pointer for jacPtr not allocated");
253 auto nb_base_functions =
255 auto nb_diff_base_functions =
259 if (nb_diff_base_functions) {
260 if (nb_diff_base_functions != nb_base_functions)
262 "Wrong number base functions %ld != %ld",
263 nb_diff_base_functions, nb_base_functions);
264 if (data.
getDiffN(base).size1() != nb_gauss_pts)
266 "Wrong number integration points");
270 if (nb_gauss_pts && nb_base_functions) {
272 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
275 double *t_transformed_n_ptr = &*
piolaN.data().begin();
278 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
283 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
284 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
285 const double a = 1. / t_det;
286 t_transformed_n(
i) =
a * (t_jac(
i,
k) * t_n(
k));
297 if (nb_gauss_pts && nb_diff_base_functions) {
302 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
304 t_transformed_diff_n(t_transformed_diff_n_ptr,
305 &t_transformed_diff_n_ptr[
HVEC0_1],
306 &t_transformed_diff_n_ptr[
HVEC0_2],
307 &t_transformed_diff_n_ptr[
HVEC1_0],
308 &t_transformed_diff_n_ptr[
HVEC1_1],
309 &t_transformed_diff_n_ptr[
HVEC1_2],
310 &t_transformed_diff_n_ptr[
HVEC2_0],
311 &t_transformed_diff_n_ptr[
HVEC2_1],
312 &t_transformed_diff_n_ptr[
HVEC2_2]);
315 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
316 for (
unsigned int bb = 0; bb != nb_diff_base_functions; ++bb) {
317 const double a = 1. / t_det;
318 t_transformed_diff_n(
i,
k) =
a * t_diff_n(
i,
k);
320 ++t_transformed_diff_n;
345 unsigned int nb_gauss_pts = data.
getN(base).size1();
346 unsigned int nb_base_functions = data.
getN(base).size2() / 3;
347 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
351 double *t_transformed_n_ptr = &*
piolaN.data().begin();
354 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
356 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
358 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[
HVEC0_1],
359 &t_transformed_diff_n_ptr[
HVEC0_2], &t_transformed_diff_n_ptr[
HVEC1_0],
360 &t_transformed_diff_n_ptr[
HVEC1_1], &t_transformed_diff_n_ptr[
HVEC1_2],
361 &t_transformed_diff_n_ptr[
HVEC2_0], &t_transformed_diff_n_ptr[
HVEC2_1],
362 &t_transformed_diff_n_ptr[
HVEC2_2]);
366 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
367 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
368 t_transformed_n(
i) = t_inv_jac(
k,
i) * t_n(
k);
369 t_transformed_diff_n(
i,
k) = t_inv_jac(
j,
i) * t_diff_n(
j,
k);
373 ++t_transformed_diff_n;
386 boost::shared_ptr<MatrixDouble> jac_ptr)
396 const auto nb_gauss_pts = getGaussPts().size2();
398 jacPtr->resize(4, nb_gauss_pts,
false);
404 auto t_t1 = getFTensor1Tangent1AtGaussPts();
405 auto t_t2 = getFTensor1Tangent2AtGaussPts();
409 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
410 t_jac(
i, N0) = t_t1(
i);
411 t_jac(
i, N1) = t_t2(
i);
429 size_t nb_gauss_pts = getGaussPts().size2();
430 jacPtr->resize(9, nb_gauss_pts,
false);
432 auto t_t1 = getFTensor1Tangent1AtGaussPts();
433 auto t_t2 = getFTensor1Tangent2AtGaussPts();
434 auto t_normal = getFTensor1NormalsAtGaussPts();
441 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
443 t_jac(
i, N0) = t_t1(
i);
444 t_jac(
i, N1) = t_t2(
i);
446 const double l = sqrt(t_normal(
j) * t_normal(
j));
448 t_jac(
i, N2) = t_normal(
i) /
l;
463 if (moab::CN::Dimension(type) != 2)
466 auto get_normals_ptr = [&]() {
473 auto apply_transform_nonlinear_geometry = [&](
auto base,
auto nb_gauss_pts,
474 auto nb_base_functions) {
477 auto ptr = get_normals_ptr();
479 &ptr[0], &ptr[1], &ptr[2]);
482 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
483 const auto l2 = t_normal(
i) * t_normal(
i);
484 for (
int bb = 0; bb != nb_base_functions; ++bb) {
485 const auto v = t_base(0);
486 t_base(
i) = (
v / l2) * t_normal(
i);
497 const auto &base_functions = data.
getN(base);
498 const auto nb_gauss_pts = base_functions.size1();
502 const auto nb_base_functions = base_functions.size2() / 3;
503 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
515 const auto type_dim = moab::CN::Dimension(type);
516 if (type_dim != 1 && type_dim != 2)
529 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
531 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
533 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
539 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
541 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
543 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
551 auto &baseN = data.
getN(base);
552 auto &diffBaseN = data.
getDiffN(base);
554 int nb_dofs = baseN.size2() / 3;
555 int nb_gauss_pts = baseN.size1();
557 piolaN.resize(baseN.size1(), baseN.size2());
558 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
560 if (nb_dofs > 0 && nb_gauss_pts > 0) {
569 t_transformed_diff_h_curl(
580 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
583 for (
int ll = 0; ll != nb_dofs; ll++) {
584 t_transformed_h_curl(
i) = t_inv_m(
j,
i) * t_h_curl(
j);
585 t_transformed_diff_h_curl(
i,
k) = t_inv_m(
j,
i) * t_diff_h_curl(
j,
k);
587 ++t_transformed_h_curl;
589 ++t_transformed_diff_h_curl;
596 if (cc != nb_gauss_pts * nb_dofs)
621 "Wrong number of integration points %ld!=%ld",
624 auto get_base_at_pts = [&](
auto base) {
631 auto get_tangent_ptr = [&]() {
639 auto get_tangent_at_pts = [&]() {
640 auto ptr = get_tangent_ptr();
642 &ptr[0], &ptr[1], &ptr[2]);
646 auto calculate_squared_edge_length = [&]() {
647 std::vector<double> l1;
649 l1.resize(nb_gauss_pts);
650 auto t_m_at_pts = get_tangent_at_pts();
651 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
652 l1[gg] = t_m_at_pts(
i) * t_m_at_pts(
i);
659 auto l1 = calculate_squared_edge_length();
664 const auto nb_dofs = data.
getN(base).size2() / 3;
666 if (nb_gauss_pts && nb_dofs) {
667 auto t_h_curl = get_base_at_pts(base);
669 auto t_m_at_pts = get_tangent_at_pts();
670 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
671 const double l0 = l1[gg];
672 for (
int ll = 0; ll != nb_dofs; ++ll) {
673 const double val = t_h_curl(0);
674 const double a = val / l0;
675 t_h_curl(
i) = t_m_at_pts(
i) *
a;
683 if (cc != nb_gauss_pts * nb_dofs)
694 boost::shared_ptr<VectorDouble> det_jac_ptr,
695 boost::shared_ptr<double> scale_det_ptr)
696 :
OP(space), fieldSpace(space), fieldBase(base), detJacPtr(det_jac_ptr),
697 scaleDetPtr(scale_det_ptr) {
708 double scale_det = 1.0;
712 auto scale = [&](
auto base) {
714 auto &base_fun = data.
getN(base);
715 if (base_fun.size2()) {
718 const auto nb_int_pts = base_fun.size1();
720 if (nb_int_pts != det_vec.size())
722 "Number of integration pts in detJacPtr does not mush "
723 "number of integration pts in base function");
725 for (
auto gg = 0; gg != nb_int_pts; ++gg) {
726 auto row = ublas::matrix_row<MatrixDouble>(base_fun, gg);
727 row *= scale_det / det_vec[gg];
730 auto &diff_base_fun = data.
getDiffN(base);
731 if (diff_base_fun.size2()) {
732 for (
auto gg = 0; gg != nb_int_pts; ++gg) {
733 auto row = ublas::matrix_row<MatrixDouble>(diff_base_fun, gg);
734 row *= scale_det / det_vec[gg];
776 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
777 std::vector<FieldSpace> spaces, std::string geom_field_name,
778 boost::shared_ptr<MatrixDouble> jac_ptr,
779 boost::shared_ptr<VectorDouble> det_ptr,
780 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
784 jac_ptr = boost::make_shared<MatrixDouble>();
786 det_ptr = boost::make_shared<VectorDouble>();
788 inv_jac_ptr = boost::make_shared<MatrixDouble>();
790 if (geom_field_name.empty()) {
805 for (
auto s : spaces) {
832 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
833 std::vector<FieldSpace> spaces, std::string geom_field_name,
834 boost::shared_ptr<MatrixDouble> jac_ptr,
835 boost::shared_ptr<VectorDouble> det_ptr,
836 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
840 jac_ptr = boost::make_shared<MatrixDouble>();
842 det_ptr = boost::make_shared<VectorDouble>();
844 inv_jac_ptr = boost::make_shared<MatrixDouble>();
846 if (geom_field_name.empty()) {
858 for (
auto s : spaces) {
888 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
889 std::vector<FieldSpace> spaces, std::string geom_field_name) {
892 if (geom_field_name.empty()) {
900 for (
auto s : spaces) {
920 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
921 std::vector<FieldSpace> spaces, std::string geom_field_name) {
924 if (geom_field_name.empty()) {
932 for (
auto s : spaces) {
947 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
948 std::vector<FieldSpace> spaces, std::string geom_field_name) {
951 if (geom_field_name.empty()) {
958 for (
auto s : spaces) {
981 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
982 std::vector<FieldSpace> spaces, std::string geom_field_name,
983 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
984 boost::shared_ptr<MatrixDouble> inv_jac) {
988 jac = boost::make_shared<MatrixDouble>();
990 det = boost::make_shared<VectorDouble>();
992 inv_jac = boost::make_shared<MatrixDouble>();
994 if (geom_field_name.empty()) {
1008 for (
auto s : spaces) {
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
static const char *const FieldSpaceNames[]
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#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.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
OpSetInvJacHcurlFaceImpl< 3 > OpSetInvJacHcurlFaceEmbeddedIn3DSpace
OpSetCovariantPiolaTransformOnFace2DImpl< 2 > OpSetCovariantPiolaTransformOnFace2D
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
OpSetContravariantPiolaTransformOnFace2DImpl< 3 > OpSetContravariantPiolaTransformOnFace2DEmbeddedIn3DSpace
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
constexpr double t
plate stiffness
FTensor::Index< 'm', 3 > m
Add operators pushing bases from local to physical configuration.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
MatrixDouble & getTangentAtGaussPts()
get tangent vector to edge curve at integration points
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, DIM > getFTensor1TangentAtGaussPts()
Get tangent vector at integration points.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives derivative)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
get hash map of derivatives base function for BB base, key is a field name
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
virtual boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
MatrixDouble & getTangent2AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
auto getFTensor1NormalsAtGaussPts()
get normal at integration points
MatrixDouble & getTangent1AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
int getFEDim() const
Get dimension of finite element.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Calculate HO coordinates at gauss points.
Calculate jacobian for face element.
boost::shared_ptr< MatrixDouble > jacPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHOJacForVolume(boost::shared_ptr< MatrixDouble > jac_ptr)
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Calculate tangent vector on edge form HO geometry approximation.
Make Hdiv space from Hcurl space in 2d.
boost::shared_ptr< VectorDouble > detJacPtr
FieldApproximationBase fieldBase
boost::shared_ptr< double > scaleDetPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpScaleBaseBySpaceInverseOfMeasure(const FieldSpace space, const FieldApproximationBase base, boost::shared_ptr< VectorDouble > det_jac_ptr, boost::shared_ptr< double > scale_det_ptr=nullptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Set inverse jacobian to base functions.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
boost::shared_ptr< MatrixDouble > invJacPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MatrixDouble diffHdivInvJac
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Modify integration weights on face to take in account higher-order geometry.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Set inverse jacobian to base functions.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > detPtr
MoFEMErrorCode applyTransform(MatrixDouble &diff_n)
Apply transformation to the input matrix.
VectorDouble & getCoords()
nodal coordinates
FTensor::Tensor2< double *, 3, 3 > & getInvJac()
get element inverse Jacobian
Volume finite element base.