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)
43 "Number of vertex coordinates and base functions is inconsistent "
45 coords.size() / 3, nb_base_functions);
48 jacPtr->resize(9, nb_gauss_pts,
false);
50 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
57 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
59 &coords[0], &coords[1], &coords[2]);
60 for (
size_t bb = 0; bb != nb_base_functions; ++bb) {
61 t_jac(
i,
j) += t_coords(
i) * (t_vol_inv_jac(
k,
j) * t_diff_base(
k));
81 return applyTransform<3, 3, 3, 3>(diff_n);
93 CHKERR transform_base(*(
m.second));
98 CHKERR transform_base(*ptr);
122 data.
getDiffN(base).size2(),
false);
124 unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
125 unsigned int nb_base_functions = data.
getDiffN(base).size2() / 9;
126 if (nb_base_functions == 0)
132 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
137 auto t_inv_jac = getFTensor2FromMat<3, 3>(*
invJacPtr);
138 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
139 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
140 t_inv_diff_n(
i,
j) = t_inv_jac(
k,
j) * t_diff_n(
i,
k);
165 for (
size_t gg = 0; gg != nb_int_pts; ++gg) {
166 t_w *= sqrt(t_normal(
i) * t_normal(
i)) /
a;
172 "Number of rows in getNormalsAtGaussPts should be equal to "
173 "number of integration points, but is not, i.e. %d != %d",
190 for (
size_t gg = 0; gg != nb_int_pts; ++gg) {
191 t_w *= sqrt(t_tangent(
i) * t_tangent(
i)) /
a;
197 "Number of rows in getTangentAtGaussPts should be equal to "
198 "number of integration points, but is not, i.e. %d != %d",
209 const auto nb_integration_pts =
detPtr->size();
214 "Inconsistent number of data points");
219 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
240 "Pointer for detPtr not allocated");
243 "Pointer for jacPtr not allocated");
251 auto nb_base_functions =
253 auto nb_diff_base_functions =
257 if (nb_diff_base_functions) {
258 if (nb_diff_base_functions != nb_base_functions)
260 "Wrong number base functions %d != %d", nb_diff_base_functions,
262 if (data.
getDiffN(base).size1() != nb_gauss_pts)
264 "Wrong number integration points");
268 if (nb_gauss_pts && nb_base_functions) {
270 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
273 double *t_transformed_n_ptr = &*
piolaN.data().begin();
276 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
279 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
281 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
282 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
283 const double a = 1. / t_det;
284 t_transformed_n(
i) =
a * t_jac(
i,
k) * t_n(
k);
295 if (nb_gauss_pts && nb_diff_base_functions) {
300 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
302 t_transformed_diff_n(t_transformed_diff_n_ptr,
303 &t_transformed_diff_n_ptr[
HVEC0_1],
304 &t_transformed_diff_n_ptr[
HVEC0_2],
305 &t_transformed_diff_n_ptr[
HVEC1_0],
306 &t_transformed_diff_n_ptr[
HVEC1_1],
307 &t_transformed_diff_n_ptr[
HVEC1_2],
308 &t_transformed_diff_n_ptr[
HVEC2_0],
309 &t_transformed_diff_n_ptr[
HVEC2_1],
310 &t_transformed_diff_n_ptr[
HVEC2_2]);
313 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
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_jac(
i,
j) * t_diff_n(
j,
k);
320 ++t_transformed_diff_n;
346 unsigned int nb_gauss_pts = data.
getN(base).size1();
347 unsigned int nb_base_functions = data.
getN(base).size2() / 3;
348 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
352 double *t_transformed_n_ptr = &*
piolaN.data().begin();
355 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
357 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
359 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[
HVEC0_1],
360 &t_transformed_diff_n_ptr[
HVEC0_2], &t_transformed_diff_n_ptr[
HVEC1_0],
361 &t_transformed_diff_n_ptr[
HVEC1_1], &t_transformed_diff_n_ptr[
HVEC1_2],
362 &t_transformed_diff_n_ptr[
HVEC2_0], &t_transformed_diff_n_ptr[
HVEC2_1],
363 &t_transformed_diff_n_ptr[
HVEC2_2]);
365 auto t_inv_jac = getFTensor2FromMat<3, 3>(*
jacInvPtr);
367 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
368 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
369 t_transformed_n(
i) = t_inv_jac(
k,
i) * t_n(
k);
370 t_transformed_diff_n(
i,
k) = t_inv_jac(
j,
i) * t_diff_n(
j,
k);
374 ++t_transformed_diff_n;
387 boost::shared_ptr<MatrixDouble> jac_ptr)
397 const auto nb_gauss_pts = getGaussPts().size2();
399 jacPtr->resize(4, nb_gauss_pts,
false);
400 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
405 auto t_t1 = getFTensor1Tangent1AtGaussPts();
406 auto t_t2 = getFTensor1Tangent2AtGaussPts();
410 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
411 t_jac(
i,
N0) = t_t1(
i);
412 t_jac(
i,
N1) = t_t2(
i);
430 size_t nb_gauss_pts = getGaussPts().size2();
431 jacPtr->resize(9, nb_gauss_pts,
false);
433 auto t_t1 = getFTensor1Tangent1AtGaussPts();
434 auto t_t2 = getFTensor1Tangent2AtGaussPts();
435 auto t_normal = getFTensor1NormalsAtGaussPts();
441 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
442 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
444 t_jac(
i,
N0) = t_t1(
i);
445 t_jac(
i,
N1) = t_t2(
i);
447 const double l = sqrt(t_normal(
j) * t_normal(
j));
449 t_jac(
i,
N2) = t_normal(
i) /
l;
464 if (moab::CN::Dimension(type) != 2)
467 auto get_normals_ptr = [&]() {
474 auto apply_transform_nonlinear_geometry = [&](
auto base,
auto nb_gauss_pts,
475 auto nb_base_functions) {
478 auto ptr = get_normals_ptr();
480 &ptr[0], &ptr[1], &ptr[2]);
483 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
484 const auto l2 = t_normal(
i) * t_normal(
i);
485 for (
int bb = 0; bb != nb_base_functions; ++bb) {
486 const auto v = t_base(0);
487 t_base(
i) = (
v / l2) * t_normal(
i);
498 const auto &base_functions = data.
getN(base);
499 const auto nb_gauss_pts = base_functions.size1();
503 const auto nb_base_functions = base_functions.size2() / 3;
504 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
516 const auto type_dim = moab::CN::Dimension(type);
517 if (type_dim != 1 && type_dim != 2)
530 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
532 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
534 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
540 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
542 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
544 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
552 auto &baseN = data.
getN(base);
553 auto &diffBaseN = data.
getDiffN(base);
555 int nb_dofs = baseN.size2() / 3;
556 int nb_gauss_pts = baseN.size1();
558 piolaN.resize(baseN.size1(), baseN.size2());
559 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
561 if (nb_dofs > 0 && nb_gauss_pts > 0) {
570 t_transformed_diff_h_curl(
581 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
584 for (
int ll = 0; ll != nb_dofs; ll++) {
585 t_transformed_h_curl(
i) = t_inv_m(
j,
i) * t_h_curl(
j);
586 t_transformed_diff_h_curl(
i,
k) = t_inv_m(
j,
i) * t_diff_h_curl(
j,
k);
588 ++t_transformed_h_curl;
590 ++t_transformed_diff_h_curl;
597 if (cc != nb_gauss_pts * nb_dofs)
622 "Wrong number of integration points %d!=%d",
625 auto get_base_at_pts = [&](
auto base) {
632 auto get_tangent_ptr = [&]() {
640 auto get_tangent_at_pts = [&]() {
641 auto ptr = get_tangent_ptr();
643 &ptr[0], &ptr[1], &ptr[2]);
647 auto calculate_squared_edge_length = [&]() {
648 std::vector<double> l1;
650 l1.resize(nb_gauss_pts);
651 auto t_m_at_pts = get_tangent_at_pts();
652 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
653 l1[gg] = t_m_at_pts(
i) * t_m_at_pts(
i);
660 auto l1 = calculate_squared_edge_length();
665 const auto nb_dofs = data.
getN(base).size2() / 3;
667 if (nb_gauss_pts && nb_dofs) {
668 auto t_h_curl = get_base_at_pts(base);
670 auto t_m_at_pts = get_tangent_at_pts();
671 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
672 const double l0 = l1[gg];
673 for (
int ll = 0; ll != nb_dofs; ++ll) {
674 const double val = t_h_curl(0);
675 const double a = val / l0;
676 t_h_curl(
i) = t_m_at_pts(
i) *
a;
684 if (cc != nb_gauss_pts * nb_dofs)
694 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
695 std::vector<FieldSpace> spaces, std::string geom_field_name,
696 boost::shared_ptr<MatrixDouble> jac_ptr,
697 boost::shared_ptr<VectorDouble> det_ptr,
698 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
702 jac_ptr = boost::make_shared<MatrixDouble>();
704 det_ptr = boost::make_shared<VectorDouble>();
706 inv_jac_ptr = boost::make_shared<MatrixDouble>();
708 if (geom_field_name.empty()) {
723 for (
auto s : spaces) {
750 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
751 std::vector<FieldSpace> spaces, std::string geom_field_name) {
754 if (geom_field_name.empty()) {
762 for (
auto s : spaces) {
782 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
783 std::vector<FieldSpace> spaces, std::string geom_field_name) {
786 if (geom_field_name.empty()) {
794 for (
auto s : spaces) {
809 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
810 std::vector<FieldSpace> spaces, std::string geom_field_name) {
813 if (geom_field_name.empty()) {
820 for (
auto s : spaces) {
841 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
842 std::vector<FieldSpace> spaces, std::string geom_field_name,
843 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
844 boost::shared_ptr<MatrixDouble> inv_jac) {
848 jac = boost::make_shared<MatrixDouble>();
850 det = boost::make_shared<VectorDouble>();
852 inv_jac = boost::make_shared<MatrixDouble>();
854 if (geom_field_name.empty()) {
868 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 .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ 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< 'm', SPACE_DIM > m
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< 2 > OpSetInvJacHcurlFace
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.
OpSetCovariantPiolaTransformOnFace2DImpl< 2 > OpSetCovariantPiolaTransformOnFace2D
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr double t
plate stiffness
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
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 direvarives 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)
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives direvatie)
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 filed 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.
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.
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
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 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.
VectorDouble & getCoords()
nodal coordinates
FTensor::Tensor2< double *, 3, 3 > & getInvJac()
get element inverse Jacobian
Volume finite element base.