95 {
96
97
99
100 auto core_log = logging::core::get();
101 core_log->add_sink(
105 core_log->add_sink(
109
110 try {
111
112
113 moab::Core moab_core;
114 moab::Interface &moab = moab_core;
115
116
119
120
121 DMType dm_name = "DMMOFEM";
123
124
127 {
128
130
132
139
143
144
145
146
147
148
150 CHKERR m_field.
get_moab().get_entities_by_type(0, MBVERTEX, verts);
155
156
158
160
163
164 auto dm = simple_interface->
getDM();
165
166 auto domain_fe =
167 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
168 auto boundary_fe =
169 boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
170
171
172 domain_fe->getRuleHook =
VolRule();
173 boundary_fe->getRuleHook =
FaceRule();
174
175 int ghosts[] = {0};
176
181
182
183 auto material_grad_mat = boost::make_shared<MatrixDouble>();
184 auto material_det_vec = boost::make_shared<VectorDouble>();
185
186 domain_fe->getOpPtrVector().push_back(
188 material_grad_mat));
189 domain_fe->getOpPtrVector().push_back(
191 domain_fe->getOpPtrVector().push_back(
193 domain_fe->getOpPtrVector().push_back(
195 domain_fe->getOpPtrVector().push_back(
196 new OpVolume(
"MESH_NODE_POSITIONS", vol));
197
198 boundary_fe->getOpPtrVector().push_back(
200 boundary_fe->getOpPtrVector().push_back(
201 new OpFace(
"MESH_NODE_POSITIONS", surf_vol));
202
204 domain_fe);
205
207 boundary_fe);
208
209 auto skeleton_fe = boost::make_shared<FEMethod>();
219 skeleton_fe->x = x;
220 skeleton_fe->x_t = x_t;
221 skeleton_fe->x_tt = x_tt;
222
223 skeleton_fe->preProcessHook = [&]() {
225 if (f != skeleton_fe->ts_F)
227 if (A != skeleton_fe->ts_A)
229 if (
B != skeleton_fe->ts_B)
231 if (x != skeleton_fe->ts_u)
233 if (x_t != skeleton_fe->ts_u_t)
235 if (x_tt != skeleton_fe->ts_u_tt)
238 };
239
240 skeleton_fe->postProcessHook = []() { return 0; };
241 skeleton_fe->operatorHook = []() { return 0; };
242
244 skeleton_fe);
245
246
247 CHKERR VecAssemblyBegin(vol);
248 CHKERR VecAssemblyEnd(vol);
249 CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
250 CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
251 CHKERR VecAssemblyBegin(surf_vol);
252 CHKERR VecAssemblyEnd(surf_vol);
253 CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
254 CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
256 double *a_vol;
257 CHKERR VecGetArray(vol, &a_vol);
258 double *a_surf_vol;
259 CHKERR VecGetArray(surf_vol, &a_surf_vol);
260 MOFEM_LOG(
"TEST", Sev::inform) <<
"Volume = " << a_vol[0];
261 MOFEM_LOG(
"TEST", Sev::inform) <<
"Surf Volume = " << a_surf_vol[0];
262 if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
264 }
265 CHKERR VecRestoreArray(vol, &a_vol);
266 CHKERR VecRestoreArray(vol, &a_surf_vol);
267 }
268
270 using ForcesAndSourcesCore::ForcesAndSourcesCore;
273 const auto type = numeredEntFiniteElementPtr->getEntType();
274 if (type != MBENTITYSET) {
276 "Expected entity set as a finite element");
277 }
280 }
283 };
284
285 auto fe_ptr = boost::make_shared<MeshsetFE>(m_field);
286
287 struct OpMeshset : ForcesAndSourcesCore::UserDataOperator {
288 using OP = ForcesAndSourcesCore::UserDataOperator;
289 using OP::OP;
290
296 }
297
298 MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
299 EntityType col_type,
306 }
307 };
308
309 fe_ptr->getOpPtrVector().push_back(
310 new OpMeshset("GLOBAL", OpMeshset::OPROW));
311 fe_ptr->getOpPtrVector().push_back(
312 new OpMeshset("GLOBAL", "GLOBAL", OpMeshset::OPROWCOL));
313
318 }
319 }
321
322
324
325 return 0;
326}
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
@ NOFIELD
scalar or vector of scalars describe (no true field)
#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.
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 createDMVector(DM dm)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createGhostVector(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_size() const =0
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
virtual MoFEMErrorCode operator()()
function is run for every finite element
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
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.
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.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode addMeshsetField(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 meshset field.
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
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 getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
std::string & getMeshsetFEName()
Get the Skeleton FE Name.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
auto & getMeshsetFiniteElementEntities()
Get the Domain Fields.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
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 getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Set integration rule to volume elements.