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",
184 const auto nb_integration_pts =
detPtr->size();
189 "Inconsistent number of data points");
194 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
215 "Pointer for detPtr not allocated");
218 "Pointer for jacPtr not allocated");
225 auto nb_gauss_pts = data.
getN(base).size1();
226 auto nb_base_functions = data.
getN(base).size2() / 3;
229 if (data.
getDiffN(base).size1() != nb_gauss_pts)
231 "Wrong number integration points");
233 if (data.
getDiffN(base).size2() / 9 != nb_base_functions)
235 "Wrong number base functions %d != %d",
236 data.
getDiffN(base).size2(), nb_base_functions);
239 if (nb_gauss_pts && nb_base_functions) {
241 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
245 double *t_transformed_n_ptr = &*
piolaN.data().begin();
248 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
250 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
252 t_transformed_diff_n(t_transformed_diff_n_ptr,
253 &t_transformed_diff_n_ptr[
HVEC0_1],
254 &t_transformed_diff_n_ptr[
HVEC0_2],
255 &t_transformed_diff_n_ptr[
HVEC1_0],
256 &t_transformed_diff_n_ptr[
HVEC1_1],
257 &t_transformed_diff_n_ptr[
HVEC1_2],
258 &t_transformed_diff_n_ptr[
HVEC2_0],
259 &t_transformed_diff_n_ptr[
HVEC2_1],
260 &t_transformed_diff_n_ptr[
HVEC2_2]);
263 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
265 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
266 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
267 const double a = 1. / t_det;
268 t_transformed_n(
i) =
a * t_jac(
i,
k) * t_n(
k);
269 t_transformed_diff_n(
i,
k) =
a * t_jac(
i,
j) * t_diff_n(
j,
k);
273 ++t_transformed_diff_n;
300 unsigned int nb_gauss_pts = data.
getN(base).size1();
301 unsigned int nb_base_functions = data.
getN(base).size2() / 3;
302 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
306 double *t_transformed_n_ptr = &*
piolaN.data().begin();
309 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
311 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
313 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[
HVEC0_1],
314 &t_transformed_diff_n_ptr[
HVEC0_2], &t_transformed_diff_n_ptr[
HVEC1_0],
315 &t_transformed_diff_n_ptr[
HVEC1_1], &t_transformed_diff_n_ptr[
HVEC1_2],
316 &t_transformed_diff_n_ptr[
HVEC2_0], &t_transformed_diff_n_ptr[
HVEC2_1],
317 &t_transformed_diff_n_ptr[
HVEC2_2]);
319 auto t_inv_jac = getFTensor2FromMat<3, 3>(*
jacInvPtr);
321 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
322 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
323 t_transformed_n(
i) = t_inv_jac(
k,
i) * t_n(
k);
324 t_transformed_diff_n(
i,
k) = t_inv_jac(
j,
i) * t_diff_n(
j,
k);
328 ++t_transformed_diff_n;
341 boost::shared_ptr<MatrixDouble> jac_ptr)
351 const auto nb_gauss_pts = getGaussPts().size2();
353 jacPtr->resize(4, nb_gauss_pts,
false);
354 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
359 auto t_t1 = getFTensor1Tangent1AtGaussPts();
360 auto t_t2 = getFTensor1Tangent2AtGaussPts();
364 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
365 t_jac(
i,
N0) = t_t1(
i);
366 t_jac(
i,
N1) = t_t2(
i);
384 size_t nb_gauss_pts = getGaussPts().size2();
385 jacPtr->resize(9, nb_gauss_pts,
false);
387 auto t_t1 = getFTensor1Tangent1AtGaussPts();
388 auto t_t2 = getFTensor1Tangent2AtGaussPts();
389 auto t_normal = getFTensor1NormalsAtGaussPts();
395 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
396 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
398 t_jac(
i,
N0) = t_t1(
i);
399 t_jac(
i,
N1) = t_t2(
i);
401 const double l = sqrt(t_normal(
j) * t_normal(
j));
403 t_jac(
i,
N2) = t_normal(
i) /
l;
418 if (moab::CN::Dimension(type) != 2)
421 auto get_normals_ptr = [&]() {
428 auto apply_transform_nonlinear_geometry = [&](
auto base,
auto nb_gauss_pts,
429 auto nb_base_functions) {
432 auto ptr = get_normals_ptr();
434 &ptr[0], &ptr[1], &ptr[2]);
437 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
438 const auto l2 = t_normal(
i) * t_normal(
i);
439 for (
int bb = 0; bb != nb_base_functions; ++bb) {
440 const auto v = t_base(0);
441 t_base(
i) = (
v / l2) * t_normal(
i);
452 const auto &base_functions = data.
getN(base);
453 const auto nb_gauss_pts = base_functions.size1();
457 const auto nb_base_functions = base_functions.size2() / 3;
458 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
470 const auto type_dim = moab::CN::Dimension(type);
471 if (type_dim != 1 && type_dim != 2)
484 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
486 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
488 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
494 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
496 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
498 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
506 auto &baseN = data.
getN(base);
507 auto &diffBaseN = data.
getDiffN(base);
509 int nb_dofs = baseN.size2() / 3;
510 int nb_gauss_pts = baseN.size1();
512 piolaN.resize(baseN.size1(), baseN.size2());
513 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
515 if (nb_dofs > 0 && nb_gauss_pts > 0) {
524 t_transformed_diff_h_curl(
535 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
538 for (
int ll = 0; ll != nb_dofs; ll++) {
539 t_transformed_h_curl(
i) = t_inv_m(
j,
i) * t_h_curl(
j);
540 t_transformed_diff_h_curl(
i,
k) = t_inv_m(
j,
i) * t_diff_h_curl(
j,
k);
542 ++t_transformed_h_curl;
544 ++t_transformed_diff_h_curl;
551 if (cc != nb_gauss_pts * nb_dofs)
576 "Wrong number of integration points %d!=%d",
579 auto get_base_at_pts = [&](
auto base) {
586 auto get_tangent_ptr = [&]() {
594 auto get_tangent_at_pts = [&]() {
595 auto ptr = get_tangent_ptr();
597 &ptr[0], &ptr[1], &ptr[2]);
601 auto calculate_squared_edge_length = [&]() {
602 std::vector<double> l1;
604 l1.resize(nb_gauss_pts);
605 auto t_m_at_pts = get_tangent_at_pts();
606 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
607 l1[gg] = t_m_at_pts(
i) * t_m_at_pts(
i);
614 auto l1 = calculate_squared_edge_length();
619 const auto nb_dofs = data.
getN(base).size2() / 3;
621 if (nb_gauss_pts && nb_dofs) {
622 auto t_h_curl = get_base_at_pts(base);
624 auto t_m_at_pts = get_tangent_at_pts();
625 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
626 const double l0 = l1[gg];
627 for (
int ll = 0; ll != nb_dofs; ++ll) {
628 const double val = t_h_curl(0);
629 const double a = val / l0;
630 t_h_curl(
i) = t_m_at_pts(
i) *
a;
638 if (cc != nb_gauss_pts * nb_dofs)
648 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
649 std::vector<FieldSpace> spaces, std::string geom_field_name,
650 boost::shared_ptr<MatrixDouble> jac_ptr,
651 boost::shared_ptr<VectorDouble> det_ptr,
652 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
656 jac_ptr = boost::make_shared<MatrixDouble>();
658 det_ptr = boost::make_shared<VectorDouble>();
660 inv_jac_ptr = boost::make_shared<MatrixDouble>();
662 if (geom_field_name.empty()) {
677 for (
auto s : spaces) {
704 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
705 std::vector<FieldSpace> spaces, std::string geom_field_name) {
708 if (geom_field_name.empty()) {
716 for (
auto s : spaces) {
736 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
737 std::vector<FieldSpace> spaces, std::string geom_field_name) {
740 if (geom_field_name.empty()) {
748 for (
auto s : spaces) {
763 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
764 std::vector<FieldSpace> spaces, std::string geom_field_name) {
767 if (geom_field_name.empty()) {
774 for (
auto s : spaces) {
795 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
796 std::vector<FieldSpace> spaces, std::string geom_field_name,
797 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
798 boost::shared_ptr<MatrixDouble> inv_jac) {
802 jac = boost::make_shared<MatrixDouble>();
804 det = boost::make_shared<VectorDouble>();
806 inv_jac = boost::make_shared<MatrixDouble>();
808 if (geom_field_name.empty()) {
822 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
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.
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
double getMeasure()
get measure of element
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.
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
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.