|
| v0.14.0
|
Go to the documentation of this file.
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)),
30 auto get_rule_by_type = [&]() {
33 return getRule(order_row + 1, order_col + 1, order_data + 1);
35 return getRule(order_row, order_col, order_data);
39 const int rule = get_rule_by_type();
41 auto calc_base_for_tet = [&]() {
43 const size_t nb_gauss_pts =
gaussPts.size2();
46 base.resize(nb_gauss_pts, 4,
false);
47 diff_base.resize(nb_gauss_pts, 12,
false);
50 double *diff_shape_ptr = &*diff_base.data().begin();
51 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
52 for (
int nn = 0; nn != 4; ++nn) {
53 for (
int dd = 0;
dd != 3; ++
dd, ++diff_shape_ptr) {
61 auto calc_base_for_hex = [&]() {
63 const size_t nb_gauss_pts =
gaussPts.size2();
66 base.resize(nb_gauss_pts, 8,
false);
67 diff_base.resize(nb_gauss_pts, 24,
false);
68 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
108 auto set_integration_pts_for_hex = [&]() {
115 auto set_integration_pts_for_tet = [&]() {
126 gaussPts.resize(4, nb_gauss_pts,
false);
127 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
129 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
131 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
140 base.resize(nb_gauss_pts, 4,
false);
141 diff_base.resize(nb_gauss_pts, 12,
false);
142 double *shape_ptr = &*base.data().begin();
143 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
145 double *diff_shape_ptr = &*diff_base.data().begin();
146 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
147 for (
int nn = 0; nn != 4; ++nn) {
148 for (
int dd = 0;
dd != 3; ++
dd, ++diff_shape_ptr) {
164 CHKERR set_integration_pts_for_tet();
167 CHKERR set_integration_pts_for_hex();
168 CHKERR calc_base_for_hex();
172 "Element type not implemented: %d",
type);
177 const size_t nb_gauss_pts =
gaussPts.size2();
181 CHKERR calc_base_for_tet();
184 CHKERR calc_base_for_hex();
188 "Element type not implemented: %d",
type);
204 auto get_tet_t_diff_n = [&]() {
210 auto get_hex_t_diff_n = [&]() {
217 auto get_t_diff_n = [&]() {
219 return get_tet_t_diff_n();
220 return get_hex_t_diff_n();
223 auto t_diff_n = get_t_diff_n();
232 tJac(
i,
j) += t_coords(
i) * t_diff_n(
j);
259 double *shape_functions_ptr = &*shape_functions.data().begin();
260 const size_t nb_base_functions = shape_functions.size2();
261 const size_t nb_gauss_pts =
gaussPts.size2();
268 shape_functions_ptr);
269 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
272 for (
int bb = 0; bb != nb_base_functions; ++bb) {
273 t_coords_at_gauss_ptr(
i) += t_coords(
i) * t_shape_functions;
277 ++t_coords_at_gauss_ptr;
290 auto dim_type = CN::Dimension(
type);
292 auto get_data_on_ents = [&](
auto lower_dim,
auto space) {
297 for (
auto dd = dim_type;
dd >= lower_dim; --
dd) {
298 int nb_ents = moab::CN::NumSubEntities(
type,
dd);
299 for (
int ii = 0; ii != nb_ents; ++ii) {
300 auto sub_ent_type = moab::CN::SubEntityType(
type,
dd, ii);
302 auto &data_on_ent = data->dataOnEntities[sub_ent_type];
305 data->spacesOnEntities[sub_ent_type].set(space);
330 if ((data.
getDiffN(base).size1() == 4) &&
331 (data.
getDiffN(base).size2() == 3)) {
332 const size_t nb_gauss_pts =
gaussPts.size2();
333 const size_t nb_base_functions = 4;
334 new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions,
false);
335 double *new_diff_n_ptr = &*new_diff_n.data().begin();
337 new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
338 double *t_diff_n_ptr = &*data.
getDiffN(base).data().begin();
339 for (
unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
341 t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
342 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
343 t_new_diff_n(
i) = t_diff_n(
i);
348 data.
getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
350 data.
getDiffN(base).swap(new_diff_n);
358 std::array<std::bitset<LASTSPACE>, 4> spaces_by_dim{
360 std::bitset<LASTSPACE>(0),
361 std::bitset<LASTSPACE>(0),
362 std::bitset<LASTSPACE>(0),
363 std::bitset<LASTSPACE>(0)
372 spaces_by_dim[1].test(
HCURL) ||
373 spaces_by_dim[2].test(
HCURL) ||
374 spaces_by_dim[3].test(
HCURL)
381 if (spaces_by_dim[2].test(
HDIV) || spaces_by_dim[3].test(
HDIV)) {
387 for (
auto t : {MBTRI, MBTET}) {
390 d.getDiffN(base) /= 6;
398 if (spaces_by_dim[3].test(
L2)) {
410 "User operator and finite element do not work together");
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Data on single entity (This is passed as argument to DataOperator::doWork)
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
Calculate base functions on tetrahedral.
MatrixInt facesNodesOrder
order of face nodes on element
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
int getMaxRowOrder() const
Get max order of approximation for field in rows.
MoFEMErrorCode getFaceNodes(EntitiesFieldData &data) const
Get nodes on faces.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
virtual MPI_Comm & get_comm() const =0
@ L2
field with C-1 continuity
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
OpSetCovariantPiolaTransform opCovariantPiolaTransform
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
int getMaxDataOrder() const
Get max order of approximation for data fields.
EntitiesFieldData & dataHdiv
OpSetInvJacHdivAndHcurl opSetInvJacHdivAndHcurl
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
Get the entity data order.
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
int getMaxColOrder() const
Get max order of approximation for field in columns.
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
MatrixInt facesNodes
nodes on finite element faces
double zeta
Viscous hardening.
Deprecated interface functions.
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.
EntitiesFieldData & dataL2
#define CHKERR
Inline error check.
virtual moab::Interface & get_moab()=0
implementation of Data Operators for Forces and Sources
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
MatrixDouble coordsAtGaussPts
coordinated at gauss points
virtual MoFEMErrorCode transformBaseFunctions()
Transform base functions based on geometric element Jacobian.
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
const EntityHandle * conn
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
OpSetInvJacH1 opSetInvJacH1
FTensor::Tensor2< double *, 3, 3 > tJac
MoFEMErrorCode operator()()
function is run for every finite element
constexpr double t
plate stiffness
Volume finite element base.
FTensor::Index< 'i', SPACE_DIM > i
static QUAD *const QUAD_3D_TABLE[]
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
structure to get information form mofem into EntitiesFieldData
#define QUAD_3D_TABLE_SIZE
EntitiesFieldData & dataH1
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
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)
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'j', 3 > j
EntitiesFieldData & dataHcurl
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ HCURL
field with continuous tangents
@ MOFEM_DATA_INCONSISTENCY
MatrixDouble gaussPts
Matrix of integration points.
OpSetContravariantPiolaTransform opContravariantPiolaTransform
FieldApproximationBase
approximation base
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Tensor2< double *, 3, 3 > tInvJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ HDIV
field with continuous normal traction
ForcesAndSourcesCore * ptrFE
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Calculate base functions on tetrahedral.
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
get sense (orientation) of entity