12static char help[] =
"...\n\n";
35 using PostProcEle::PostProcEle;
92 PetscBool load_file = PETSC_FALSE;
96 if (load_file == PETSC_FALSE) {
103 double tet_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
105 for (
int nn = 0; nn < 4; nn++) {
106 CHKERR moab.create_vertex(&tet_coords[3 * nn], nodes[nn]);
109 CHKERR moab.create_element(MBTET, nodes, 4, tet);
111 for (
auto d : {1, 2})
112 CHKERR moab.get_adjacencies(&tet, 1, d,
true, adj);
118 double tri_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
120 for (
int nn = 0; nn < 3; nn++) {
121 CHKERR moab.create_vertex(&tri_coords[3 * nn], nodes[nn]);
124 CHKERR moab.create_element(MBTRI, nodes, 3, tri);
126 CHKERR moab.get_adjacencies(&tri, 1, 1,
true, adj);
154 enum bases { AINSWORTH, AINSWORTH_LOBATTO, DEMKOWICZ, BERNSTEIN, LASBASETOP };
155 const char *list_bases[] = {
"ainsworth",
"ainsworth_lobatto",
"demkowicz",
159 PetscInt choice_base_value = AINSWORTH;
161 &choice_base_value, &flg);
162 if (flg != PETSC_TRUE)
165 if (choice_base_value == AINSWORTH)
167 if (choice_base_value == AINSWORTH_LOBATTO)
169 else if (choice_base_value == DEMKOWICZ)
171 else if (choice_base_value == BERNSTEIN)
174 enum spaces { H1SPACE, L2SPACE, HCURLSPACE, HDIVSPACE, LASBASETSPACE };
175 const char *list_spaces[] = {
"h1",
"l2",
"hcurl",
"hdiv"};
176 PetscInt choice_space_value = H1SPACE;
178 LASBASETSPACE, &choice_space_value, &flg);
179 if (flg != PETSC_TRUE)
182 if (choice_space_value == H1SPACE)
184 else if (choice_space_value == L2SPACE)
186 else if (choice_space_value == HCURLSPACE)
188 else if (choice_space_value == HDIVSPACE)
231 auto post_proc_fe = boost::make_shared<MyPostProc>(
mField);
232 post_proc_fe->generateReferenceElementMesh();
237 auto jac_ptr = boost::make_shared<MatrixDouble>();
238 post_proc_fe->getOpPtrVector().push_back(
241 post_proc_fe->getOpPtrVector().push_back(
254 auto u_ptr = boost::make_shared<VectorDouble>();
255 post_proc_fe->getOpPtrVector().push_back(
257 post_proc_fe->getOpPtrVector().push_back(
261 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
281 auto u_ptr = boost::make_shared<MatrixDouble>();
282 post_proc_fe->getOpPtrVector().push_back(
285 post_proc_fe->getOpPtrVector().push_back(
289 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
307 auto scale_tag_val = [&]() {
309 auto &post_proc_mesh = post_proc_fe->getPostProcMesh();
311 CHKERR post_proc_mesh.get_entities_by_type(0, MBVERTEX, nodes);
313 CHKERR post_proc_mesh.tag_get_handle(
"U",
th);
315 CHKERR post_proc_mesh.tag_get_length(
th, length);
316 std::vector<double> data(nodes.size() * length);
317 CHKERR post_proc_mesh.tag_get_data(
th, nodes, &*data.begin());
319 for (
int i = 0;
i != nodes.size(); ++
i) {
321 for (
int d = 0;
d != length; ++
d)
322 v += pow(data[length *
i + d], 2);
324 max_v = std::max(max_v,
v);
328 CHKERR post_proc_mesh.tag_set_data(
th, nodes, &*data.begin());
333 auto dofs_ptr = mField.get_dofs();
335 for (
auto dof_ptr : (*dofs_ptr)) {
336 MOFEM_LOG(
"PLOTBASE", Sev::verbose) << *dof_ptr;
337 auto &val =
const_cast<double &
>(dof_ptr->getFieldData());
341 CHKERR post_proc_fe->writeFile(
342 "out_base_dof_" + boost::lexical_cast<std::string>(nb) +
".h5m");
343 CHKERR post_proc_fe->getPostProcMesh().delete_mesh();
356int main(
int argc,
char *argv[]) {
364 DMType dm_name =
"DMMOFEM";
369 auto core_log = logging::core::get();
376 moab::Core mb_instance;
377 moab::Interface &moab = mb_instance;
400 moab::Interface &moab_ref = core_ref;
402 char ref_mesh_file_name[255];
405 strcpy(ref_mesh_file_name,
"ref_mesh2d.h5m");
407 strcpy(ref_mesh_file_name,
"ref_mesh3d.h5m");
410 "Dimension not implemented");
414 CHKERR moab_ref.load_file(ref_mesh_file_name, 0,
"");
422 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
423 CHKERR moab_ref.add_entities(meshset, elems);
424 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
425 CHKERR moab_ref.delete_entities(&meshset, 1);
429 CHKERR moab_ref.get_connectivity(elems, elem_nodes,
false);
432 std::map<EntityHandle, int> nodes_pts_map;
435 gaussPts.resize(
SPACE_DIM + 1, elem_nodes.size(),
false);
437 Range::iterator nit = elem_nodes.begin();
438 for (
int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
440 CHKERR moab_ref.get_coords(&*nit, 1, coords);
441 for (
auto d : {0, 1, 2})
442 gaussPts(d, gg) = coords[d];
443 nodes_pts_map[*nit] = gg;
455 Range::iterator tit = elems.begin();
456 for (
int tt = 0; tit != elems.end(); ++tit, ++tt) {
459 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
460 for (
int nn = 0; nn != num_nodes; ++nn) {
461 refEleMap(tt, nn) = nodes_pts_map[conn[nn]];
471 const int num_nodes = gaussPts.size2();
475 switch (numeredEntFiniteElementPtr->getEntType()) {
479 &gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
483 for (
int gg = 0; gg != num_nodes; gg++) {
484 double ksi = gaussPts(0, gg);
485 double eta = gaussPts(1, gg);
495 &gaussPts(0, 0), &gaussPts(1, 0),
496 &gaussPts(2, 0), num_nodes);
500 for (
int gg = 0; gg != num_nodes; gg++) {
501 double ksi = gaussPts(0, gg);
502 double eta = gaussPts(1, gg);
503 double zeta = gaussPts(2, gg);
516 "Not implemented element type");
522 ReadUtilIface *iface;
523 CHKERR getPostProcMesh().query_interface(iface);
525 std::vector<double *> arrays;
530 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
532 mapGaussPts.resize(gaussPts.size2());
533 for (
int gg = 0; gg != num_nodes; ++gg)
534 mapGaussPts[gg] = startv + gg;
537 int def_in_the_loop = -1;
538 CHKERR getPostProcMesh().tag_get_handle(
"NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
539 th, MB_TAG_CREAT | MB_TAG_SPARSE,
545 const int num_nodes_on_ele =
refEleMap.size2();
552 CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
555 CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
559 "Dimension not implemented");
563 for (
unsigned int tt = 0; tt !=
refEleMap.size1(); ++tt) {
564 for (
int nn = 0; nn != num_nodes_on_ele; ++nn)
565 conn[num_nodes_on_ele * tt + nn] = mapGaussPts[
refEleMap(tt, nn)];
570 CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
572 auto physical_elements =
Range(starte, starte + num_el - 1);
573 CHKERR getPostProcMesh().tag_clear_data(
th, physical_elements, &(nInTheLoop));
575 EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
579 mField.get_moab().get_connectivity(fe_ent, conn, fe_num_nodes,
true);
580 coords.resize(3 * fe_num_nodes,
false);
581 CHKERR mField.get_moab().get_coords(conn, fe_num_nodes, &coords[0]);
590 arrays[0], arrays[1], arrays[2]);
591 const double *t_coords_ele_x = &coords[0];
592 const double *t_coords_ele_y = &coords[1];
593 const double *t_coords_ele_z = &coords[2];
594 for (
int gg = 0; gg != num_nodes; ++gg) {
596 t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
598 for (
int nn = 0; nn != fe_num_nodes; ++nn) {
599 t_coords(
i) += t_n * t_ele_coords(
i);
600 for (
auto ii : {0, 1, 2})
601 if (std::abs(t_coords(ii)) < std::numeric_limits<float>::epsilon())
614 ParallelComm *pcomm_post_proc_mesh =
616 if (pcomm_post_proc_mesh != NULL)
617 delete pcomm_post_proc_mesh;
624 auto resolve_shared_ents = [&]() {
627 ParallelComm *pcomm_post_proc_mesh =
628 ParallelComm::get_pcomm(&(getPostProcMesh()),
MYPCOMM_INDEX);
629 if (pcomm_post_proc_mesh == NULL) {
632 pcomm_post_proc_mesh =
new ParallelComm(
633 &(getPostProcMesh()),
637 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
642 CHKERR resolve_shared_ents();
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#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 ...
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
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.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, 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
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
double zeta
Viscous hardening.
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode setIntegrationRules()
[Set up problem]
FieldApproximationBase base
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
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)
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.
Get vector field for H-div approximation.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
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.
void setDim(int dim)
Set the problem dimension.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
MoFEMErrorCode postProcess()
MoFEMErrorCode generateReferenceElementMesh()
ublas::matrix< int > refEleMap
MoFEMErrorCode setGaussPts(int order)
MoFEMErrorCode preProcess()
MatrixDouble shapeFunctions