v0.14.0
Loading...
Searching...
No Matches
projection_from_10node_op_volume.cpp
Go to the documentation of this file.
1
2
3#include <MoFEM.hpp>
4
5using namespace MoFEM;
6
7static char help[] = "...\n\n";
8
11
13 OpVolumeCalculation(Vec volume_vec)
15 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
16 volumeVec(volume_vec) {}
17
18 MoFEMErrorCode doWork(int row_side, EntityType row_type,
21
22 if (row_type != MBVERTEX)
24 double vol = 0;
26
27 int nb_gauss_pts = row_data.getN().size1();
28 for (int gg = 0; gg != nb_gauss_pts; gg++) {
29 vol += getVolume() * t_w;
30 ++t_w;
31 }
32
33 CHKERR VecSetValue(volumeVec, 0, vol, ADD_VALUES);
34
36 }
37};
38
39int main(int argc, char *argv[]) {
40
41 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
42
43 try {
44
45 moab::Core mb_instance;
46 moab::Interface &moab = mb_instance;
47 int rank;
48 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
49
50 PetscBool flg = PETSC_TRUE;
51 char mesh_file_name[255];
52#if PETSC_VERSION_GE(3, 6, 4)
53 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
54 255, &flg);
55#else
56 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
57 mesh_file_name, 255, &flg);
58#endif
59 if (flg != PETSC_TRUE) {
60 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
61 }
62
63 const char *option;
64 option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
65 CHKERR moab.load_file(mesh_file_name, 0, option);
66
67 // ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
68 // auto moab_comm_wrap =
69 // boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
70 // if (pcomm == NULL)
71 // pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
72
73 MoFEM::Core core(moab);
74 MoFEM::Interface &m_field = core;
75
76 // add filds
77 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
78 3);
79
80 // add finite elements
81 CHKERR m_field.add_finite_element("TET_ELEM");
83 "MESH_NODE_POSITIONS");
85 "MESH_NODE_POSITIONS");
87 "MESH_NODE_POSITIONS");
88
89 BitRefLevel bit_level0;
90 bit_level0.set(0);
91 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
92 0, 3, bit_level0);
93 Range tets;
94 CHKERR moab.get_entities_by_type(0, MBTET, tets, true);
95 Range edges;
96 CHKERR moab.get_entities_by_type(0, MBEDGE, edges, true);
97 CHKERR m_field.getInterface<BitRefManager>()->setElementsBitRefLevel(edges);
98
99 // add ents to field and set app. order
100 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
101 CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
102 CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
103 CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
104 CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
105
106 // add finite elements entities
107 CHKERR m_field.add_ents_to_finite_element_by_type(tets, MBTET, "TET_ELEM");
108
109 // build fields
110 CHKERR m_field.build_fields();
111 // build finite elements
113 // build adjacencies
114 CHKERR m_field.build_adjacencies(bit_level0);
115
116 Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
117 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
118
119 CHKERR DMRegister_MoFEM("DMMOFEM");
120 auto dM = createDM(m_field.get_comm(), "DMMOFEM");
121 CHKERR DMMoFEMCreateMoFEM(dM, &m_field, "TET_PROBLEM", bit_level0);
122 CHKERR DMMoFEMAddElement(dM, "TET_ELEM");
123 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_FALSE);
125
126 auto vol_vec = createVectorMPI(m_field.get_comm(),
127 m_field.get_comm_rank() ? 0 : 1, 1);
128
129 auto fe_ptr =
130 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
131 auto material_grad_mat = boost::make_shared<MatrixDouble>();
132 auto material_det_vec = boost::make_shared<VectorDouble>();
133 fe_ptr->meshPositionsFieldName = "none";
134 fe_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
135 "MESH_NODE_POSITIONS", material_grad_mat));
136 fe_ptr->getOpPtrVector().push_back(
137 new OpInvertMatrix<3>(material_grad_mat, material_det_vec, nullptr));
138 fe_ptr->getOpPtrVector().push_back(new OpSetHOWeights(material_det_vec));
139 fe_ptr->getOpPtrVector().push_back(new OpVolumeCalculation(vol_vec));
140
141 CHKERR DMoFEMLoopFiniteElements(dM, "TET_ELEM", fe_ptr);
142
143 double vol;
144 CHKERR VecAssemblyBegin(vol_vec);
145 CHKERR VecAssemblyEnd(vol_vec);
146 CHKERR VecSum(vol_vec, &vol);
147 CHKERR PetscPrintf(PETSC_COMM_WORLD, "vol = %3.8e\n", vol);
148 if (fabs(vol - 1.0) > 1e-4)
149 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Failed to pass test");
150 }
152
154}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
Definition: adol-c_atom.cpp:46
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
PetscErrorCode DMSetUp_MoFEM(DM dm)
Definition: DMMoFEM.cpp:1267
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
char mesh_file_name[255]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
static char help[]
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Set inverse jacobian to base functions.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.