12 meshPositionsFieldName(
"MESH_NODE_POSITIONS"), coords(24), jAc(3, 3),
14 opContravariantPiolaTransform(elementMeasure, jAc),
15 opCovariantPiolaTransform(
invJac), opSetInvJacHdivAndHcurl(
invJac),
16 tJac(&jAc(0, 0), &jAc(0, 1), &jAc(0, 2), &jAc(1, 0), &jAc(1, 1),
17 &jAc(1, 2), &jAc(2, 0), &jAc(2, 1), &jAc(2, 2)),
23 Interface &m_field,
const EntityType type)
34 auto get_rule_by_type = [&]() {
37 return getRule(order_row + 1, order_col + 1, order_data + 1);
39 return getRule(order_row, order_col, order_data);
43 const int rule = get_rule_by_type();
45 auto calc_base_for_tet = [&]() {
47 const size_t nb_gauss_pts =
gaussPts.size2();
50 base.resize(nb_gauss_pts, 4,
false);
51 diff_base.resize(nb_gauss_pts, 12,
false);
54 double *diff_shape_ptr = &*diff_base.data().begin();
55 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
56 for (
int nn = 0; nn != 4; ++nn) {
57 for (
int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
65 auto calc_base_for_hex = [&]() {
67 const size_t nb_gauss_pts =
gaussPts.size2();
70 base.resize(nb_gauss_pts, 8,
false);
71 diff_base.resize(nb_gauss_pts, 24,
false);
72 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
112 auto set_integration_pts_for_hex = [&]() {
119 auto set_integration_pts_for_tet = [&]() {
130 gaussPts.resize(4, nb_gauss_pts,
false);
131 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
133 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
135 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
144 base.resize(nb_gauss_pts, 4,
false);
145 diff_base.resize(nb_gauss_pts, 12,
false);
146 double *shape_ptr = &*base.data().begin();
147 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
149 double *diff_shape_ptr = &*diff_base.data().begin();
150 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
151 for (
int nn = 0; nn != 4; ++nn) {
152 for (
int dd = 0;
dd != 3; ++
dd, ++diff_shape_ptr) {
168 CHKERR set_integration_pts_for_tet();
171 CHKERR set_integration_pts_for_hex();
172 CHKERR calc_base_for_hex();
176 "Element type not implemented: %d", type);
181 const size_t nb_gauss_pts =
gaussPts.size2();
185 CHKERR calc_base_for_tet();
188 CHKERR calc_base_for_hex();
192 "Element type not implemented: %d", type);
208 auto get_tet_t_diff_n = [&]() {
214 auto get_hex_t_diff_n = [&]() {
221 auto get_t_diff_n = [&]() {
223 return get_tet_t_diff_n();
224 return get_hex_t_diff_n();
227 auto t_diff_n = get_t_diff_n();
236 tJac(
i,
j) += t_coords(
i) * t_diff_n(
j);
263 double *shape_functions_ptr = &*shape_functions.data().begin();
264 const size_t nb_base_functions = shape_functions.size2();
265 const size_t nb_gauss_pts =
gaussPts.size2();
272 shape_functions_ptr);
273 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
276 for (
int bb = 0; bb != nb_base_functions; ++bb) {
277 t_coords_at_gauss_ptr(
i) += t_coords(
i) * t_shape_functions;
281 ++t_coords_at_gauss_ptr;
294 auto dim_type = CN::Dimension(type);
296 auto get_data_on_ents = [&](
auto lower_dim,
auto space) {
301 for (
auto dd = dim_type; dd >= lower_dim; --dd) {
302 int nb_ents = moab::CN::NumSubEntities(type, dd);
303 for (
int ii = 0; ii != nb_ents; ++ii) {
304 auto sub_ent_type = moab::CN::SubEntityType(type, dd, ii);
306 auto &data_on_ent = data->dataOnEntities[sub_ent_type];
309 data->spacesOnEntities[sub_ent_type].set(space);
334 if ((data.
getDiffN(base).size1() == 4) &&
335 (data.
getDiffN(base).size2() == 3)) {
336 const size_t nb_gauss_pts =
gaussPts.size2();
337 const size_t nb_base_functions = 4;
338 new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions,
false);
339 double *new_diff_n_ptr = &*new_diff_n.data().begin();
341 new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
342 double *t_diff_n_ptr = &*data.
getDiffN(base).data().begin();
343 for (
unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
345 t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
346 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
347 t_new_diff_n(
i) = t_diff_n(
i);
352 data.
getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
354 data.
getDiffN(base).swap(new_diff_n);
362 std::array<std::bitset<LASTSPACE>, 4> spaces_by_dim{
364 std::bitset<LASTSPACE>(0),
365 std::bitset<LASTSPACE>(0),
366 std::bitset<LASTSPACE>(0),
367 std::bitset<LASTSPACE>(0)
370 for (
auto type = MBEDGE; type != MBENTITYSET; ++type) {
376 spaces_by_dim[1].test(
HCURL) ||
377 spaces_by_dim[2].test(
HCURL) ||
378 spaces_by_dim[3].test(
HCURL)
385 if (spaces_by_dim[2].test(
HDIV) || spaces_by_dim[3].test(
HDIV)) {
391 for (
auto t : {MBTRI, MBTET}) {
394 d.getDiffN(base) /= 6;
402 if (spaces_by_dim[3].test(
L2)) {
414 "User operator and finite element do not work together");
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#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
#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 MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
FTensor::Index< 'j', 3 > j
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
double zeta
Viscous hardening.
constexpr double t
plate stiffness
#define QUAD_3D_TABLE_SIZE
static QUAD *const QUAD_3D_TABLE[]
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MatrixInt facesNodesOrder
order of face nodes on element
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
ForcesAndSourcesCore * ptrFE
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
int getMaxRowOrder() const
Get max order of approximation for field in rows.
EntitiesFieldData & dataHdiv
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getFaceNodes(EntitiesFieldData &data) const
Get nodes on faces.
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
EntitiesFieldData & dataHcurl
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
EntitiesFieldData & dataH1
MatrixDouble coordsAtGaussPts
coordinated at gauss points
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
get sense (orientation) of entity
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
int getMaxColOrder() const
Get max order of approximation for field in columns.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
Get the entity data order.
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
EntitiesFieldData & dataL2
Calculate base functions on tetrahedral.
Calculate base functions on tetrahedral.
Volume finite element base.
OpSetInvJacHdivAndHcurl opSetInvJacHdivAndHcurl
FTensor::Tensor2< double *, 3, 3 > tJac
VolumeElementForcesAndSourcesCore(Interface &m_field)
virtual MoFEMErrorCode transformBaseFunctions()
Transform base functions based on geometric element Jacobian.
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
MoFEMErrorCode operator()()
function is run for every finite element
OpSetCovariantPiolaTransform opCovariantPiolaTransform
FTensor::Tensor2< double *, 3, 3 > tInvJac
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
OpSetContravariantPiolaTransform opContravariantPiolaTransform
const EntityHandle * conn
OpSetInvJacH1 opSetInvJacH1