95 {
96
97
99
100 try {
101
102
103 moab::Core moab_core;
104 moab::Interface &moab = moab_core;
105
106
109
110
111 DMType dm_name = "DMMOFEM";
113
114
117 {
118
120
122
129
130
132
134
137
138 auto dm = simple_interface->
getDM();
139
140 auto domain_fe =
141 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
142 auto boundary_fe =
143 boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
144
145
146 domain_fe->getRuleHook =
VolRule();
147 boundary_fe->getRuleHook =
FaceRule();
148
149 int ghosts[] = {0};
150
155
156
157 auto material_grad_mat = boost::make_shared<MatrixDouble>();
158 auto material_det_vec = boost::make_shared<VectorDouble>();
159
160 domain_fe->meshPositionsFieldName = "none";
161 domain_fe->getOpPtrVector().push_back(
163 material_grad_mat));
165 material_grad_mat, material_det_vec, nullptr));
166 domain_fe->getOpPtrVector().push_back(
168 domain_fe->getOpPtrVector().push_back(
170 domain_fe->getOpPtrVector().push_back(
171 new OpVolume(
"MESH_NODE_POSITIONS", vol));
172
173 boundary_fe->getOpPtrVector().push_back(
175 boundary_fe->getOpPtrVector().push_back(
176 new OpFace(
"MESH_NODE_POSITIONS", surf_vol));
177
179 domain_fe);
180
182 boundary_fe);
183
184 auto skeleton_fe = boost::make_shared<FEMethod>();
194 skeleton_fe->x = x;
195 skeleton_fe->x_t = x_t;
196 skeleton_fe->x_tt = x_tt;
197
198 skeleton_fe->preProcessHook = [&]() {
200 if (f != skeleton_fe->ts_F)
202 if (A != skeleton_fe->ts_A)
204 if (
B != skeleton_fe->ts_B)
206 if (x != skeleton_fe->ts_u)
208 if (x_t != skeleton_fe->ts_u_t)
210 if (x_tt != skeleton_fe->ts_u_tt)
213 };
214
215 skeleton_fe->postProcessHook = []() { return 0; };
216 skeleton_fe->operatorHook = []() { return 0; };
217
219 skeleton_fe);
220
221
222 CHKERR VecAssemblyBegin(vol);
223 CHKERR VecAssemblyEnd(vol);
224 CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
225 CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
226 CHKERR VecAssemblyBegin(surf_vol);
227 CHKERR VecAssemblyEnd(surf_vol);
228 CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
229 CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
231 double *a_vol;
232 CHKERR VecGetArray(vol, &a_vol);
233 double *a_surf_vol;
234 CHKERR VecGetArray(surf_vol, &a_surf_vol);
235 cout << "Volume = " << a_vol[0] << endl;
236 cout << "Surf Volume = " << a_surf_vol[0] << endl;
237 if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
239 }
240 CHKERR VecRestoreArray(vol, &a_vol);
241 CHKERR VecRestoreArray(vol, &a_surf_vol);
242 }
243
244 }
245 }
247
248
250
251 return 0;
252}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
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.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
auto createSmartGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
Set integration rule to boundary elements.
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Calculate HO coordinates at gauss points.
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.
Simple interface for fast problem set-up.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
MoFEMErrorCode addDomainField(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_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode addBoundaryField(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_ZERO, int verb=-1)
Add field on boundary.
MoFEMErrorCode addSkeletonField(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_ZERO, int verb=-1)
Add field on skeleton.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.