24#include <boost/numeric/ublas/symmetric.hpp>
32#include <adolc/adolc.h>
43#if PETSC_VERSION_GE(3, 6, 0)
44#include <petsc/private/dmimpl.h>
45#include <petsc/private/pcmgimpl.h>
47#include <petsc-private/dmimpl.h>
48#include <petsc-private/pcmgimpl.h>
63 PetscFunctionReturn(0);
70 PetscFunctionReturn(0);
75 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
76 "MOFEM Shell Multi-Grid (Orders) pre-conditioner",
79 CHKERR PetscOptionsInt(
"-mofem_mg_verbose",
"nb levels of multi-grid solver",
80 "", 0, &verboseLevel, PETSC_NULL);
82 ierr = PetscOptionsEnd();
85 PetscFunctionReturn(0);
92 verb = verb > verboseLevel ? verb : verboseLevel;
95 CHKERR PetscObjectGetComm((PetscObject)
dM, &comm);
98 PetscPrintf(comm,
"set MG levels %u\n", nbLevels);
101 vector<IS> is_vec(nbLevels + 1);
102 vector<int> is_glob_size(nbLevels + 1), is_loc_size(nbLevels + 1);
104 for (
int kk = 0; kk < nbLevels; kk++) {
107 CHKERR createIsAtLevel(kk, &is_vec[kk]);
108 CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
109 CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
112 PetscSynchronizedPrintf(comm,
113 "Nb. dofs at level [ %d ] global %u local %d\n",
114 kk, is_glob_size[kk], is_loc_size[kk]);
118 if (is_glob_size[kk] == 0) {
130#if PETSC_VERSION_GE(3, 8, 0)
131 CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_PMAT);
132#warning "This is not working downgrade petsc-3.7.* or optimally petsc-3.6.*"
134 CHKERR PCMGSetGalerkin(pc, PETSC_TRUE);
142 for (
unsigned int kk = 0; kk < is_vec.size(); kk++) {
143 CHKERR destroyIsAtLevel(kk, &is_vec[kk]);
147 PetscSynchronizedFlush(comm, PETSC_STDOUT);
150 PetscFunctionReturn(0);
155 PetscFunctionReturn(0);
160 PetscFunctionReturn(0);
170 for (_IT_FENUMEREDDOF_BY_NAME_ROW_FOR_LOOP_(numeredEntFiniteElementPtr,
171 "DISPLACEMENT", dof)) {
172 if (dof->get()->getEntType() == MBVERTEX)
174 if (ent != dof->get()->getEnt()) {
175 ent = dof->get()->getEnt();
179 vector<int> adj_orders(adj_prisms.size());
183 vector<int>::iterator max_it =
184 max_element(adj_orders.begin(), adj_orders.end());
186 if (
order > max_order) {
189 "order of element can not be bigger than order of entity %d > %d",
194 if (
order < max_order) {
201 PetscFunctionReturn(0);
209 "ENERGY_OR_ENERGY_ERROR", 1, MB_TYPE_DOUBLE, thError,
210 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
213 PetscFunctionReturn(0);
218 PetscFunctionReturn(0);
223 EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
228 if (
order == maxOrder)
229 PetscFunctionReturn(0);
235 if (elem_error < eRror && order > 2) {
245 if (
order < maxOrder) {
253 PetscFunctionReturn(0);
257 int max_order,
int verb) {
261 thOrder, MB_TAG_CREAT | MB_TAG_SPARSE,
275 ParallelComm *pcomm =
278 rval = pcomm->exchange_tags(
th, prisms);
282 tags.push_back(pcomm->part_tag());
286 tags.push_back(th_global_id);
290 o1 <<
"test_" << rr <<
".vtk";
291 mField.
get_moab().write_file(o1.str().c_str(),
"VTK",
"", &meshset, 1,
292 &tags[0], tags.size());
299 PetscFunctionReturn(0);
309 for (_IT_NUMEREDDOF_ROW_BY_OWNPROC_FOR_LOOP_(problem_ptr,
311 if (dof->get()->getName() !=
"DISPLACEMENT") {
312 vecOrderDofs.push_back(dof->get()->getPetscGlobalDofIdx());
313 }
else if (dof->get()->getEntType() == MBVERTEX) {
314 vecOrderDofs.push_back(dof->get()->getPetscGlobalDofIdx());
316 if (ent != dof->get()->getEnt()) {
317 ent = dof->get()->getEnt();
321 if (dof->get()->getDofOrder() <= ent_order) {
322 vecOrderDofs.push_back(dof->get()->getPetscGlobalDofIdx());
328 PetscFunctionReturn(0);
333 PetscPrintf(PETSC_COMM_WORLD,
"SetUp sOlver ... ");
339 PetscBool same = PETSC_FALSE;
341 PetscErrorCode (*org_pc_setup_mg)(PC pc) =
pC->ops->setup;
343 PetscObjectTypeCompare((PetscObject)
pC, PCMG, &same);
345 PetscPrintf(PETSC_COMM_WORLD,
"multigrid ..\n");
356 "Expected MG pre-conditioner");
360 PetscPrintf(PETSC_COMM_WORLD,
"Solve problem ..\n");
363 PetscErrorCode (*org_reset)(PC pc) =
pC->ops->reset;
366 pC->ops->reset = org_reset;
367 pC->ops->setup = org_pc_setup_mg;
369 PetscFunctionReturn(0);
375 PetscPrintf(PETSC_COMM_WORLD,
"Get level sub matrix ...\n");
378 MAT_INITIAL_MATRIX, &sub_Aij);
391 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
392 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
393 CHKERR MatDestroy(&sub_Aij);
394 PetscFunctionReturn(0);
405 for (
int ii = 0; ii < loc_size; ii++) {
409 int procent_size = ceil(size * procent);
410 procent_size = procent_size > loc_size ? loc_size : procent_size;
413 PETSC_DETERMINE,
a0, &short_sorted);
417 CHKERR VecScatterCreateToZero(short_sorted, &ctx, &all_in_zero);
418 CHKERR VecScatterBegin(ctx, short_sorted, all_in_zero, INSERT_VALUES,
420 CHKERR VecScatterEnd(ctx, short_sorted, all_in_zero, INSERT_VALUES,
424 CHKERR VecGetArray(all_in_zero, &
a1);
425 int all_in_zero_size;
426 CHKERR VecGetSize(all_in_zero, &all_in_zero_size);
427 CHKERR PetscSortReal(all_in_zero_size,
a1);
428 int pos = (
unsigned int)ceil(size * procent);
429 pos = pos > all_in_zero_size ? all_in_zero_size - 1 : pos;
432 CHKERR VecRestoreArray(all_in_zero, &
a1);
434 CHKERR VecDestroy(&all_in_zero);
435 CHKERR VecScatterDestroy(&ctx);
436 CHKERR VecDestroy(&short_sorted);
443 CHKERR VecGhostUpdateBegin(exchange_value, INSERT_VALUES, SCATTER_FORWARD);
444 CHKERR VecGhostUpdateEnd(exchange_value, INSERT_VALUES, SCATTER_FORWARD);
445 CHKERR VecDestroy(&exchange_value);
446 for (
int ii = 0; ii < loc_size; ii++) {
453 CHKERR VecDestroy(&short_sorted);
455 PetscFunctionReturn(0);
459 int nb_ref_cycles,
int max_order,
460 double procent,
int verb) {
463 const int max_levels = 40;
464 isLevels.resize(max_levels, PETSC_NULL);
473 CHKERR PetscLayoutGetLocalSize(layout, &size);
474 CHKERR VecCreateMPI(PETSC_COMM_WORLD, size, PETSC_DETERMINE,
476 CHKERR PetscLayoutDestroy(&layout);
496 PetscPrintf(PETSC_COMM_WORLD,
"Level %d Energy = %12.9e\n", 0, energy);
509 double max_error = 0;
515 CHKERR VecGhostUpdateBegin(
E, INSERT_VALUES, SCATTER_FORWARD);
516 CHKERR VecGhostUpdateEnd(
E, INSERT_VALUES, SCATTER_FORWARD);
517 CHKERR VecNorm(
E, NORM_2, &vnorm);
524 PetscPrintf(PETSC_COMM_WORLD,
525 "Refinement %d Size = %d vec norm = %8.6e error at precent "
526 "(default 33) = %12.9e ( max error %12.9e )\n",
532 CHKERR VecNorm(
D, NORM_2, &vnorm);
537 PetscPrintf(PETSC_COMM_WORLD,
538 "Level %d size = %d vec norm = %6.4e energy = %12.9e\n",
539 rr + 1, is_size, vnorm, energy);
546 if (rr == nb_ref_cycles) {
562 PetscPrintf(PETSC_COMM_WORLD,
563 "Last level size = %d, previous level size = %d\n", s1, s2);
577 CHKERR VecCopy(sub_b, sub_d);
586 VecGetSubVector(
D,
mgLevels.back(), &sub_d);
588 CHKERR VecCopy(sub_d, sub_b);
589 VecRestoreSubVector(
D,
mgLevels.back(), &sub_d);
598 for (
int ll = 0; ll < max_levels; ll++) {
614 PetscFunctionReturn(0);
653 PetscFunctionReturn(0);
661 "PARALLEL=WRITE_PART");
663 PetscFunctionReturn(0);
674 mField.
get_moab().get_entities_by_type(meshset, MBVERTEX, nodes,
true);
678 moab::Interface::UNION);
682 moab::Interface::UNION);
684 Range tris = nodes_faces.subset_by_type(MBTRI);
690 moab::Interface::UNION);
693 quads = faces.subset_by_type(MBQUAD);
700 nodes.merge(subtract(intersect(edges_nodes, quad_nodes), tris_nodes));
701 Range::iterator nit = nodes.begin();
702 for (; nit != nodes.end(); nit++) {
703 DofEntityByNameAndEnt::iterator dit, hi_dit;
705 boost::make_tuple(
"DISPLACEMENT", *nit));
707 boost::make_tuple(
"DISPLACEMENT", *nit));
708 for (; dit != hi_dit; dit++) {
713 PetscSynchronizedPrintf(
715 "DISPLACEMENT [ %d ] %6.4e coord [ %3.4f %3.4f %3.4f ] "
716 "EntityHandle %ld rank %d\n",
717 dit->get()->getDofCoeffIdx(), dit->get()->getFieldData(), coords[0],
722 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
723 PetscFunctionReturn(0);
MoFEMErrorCode DMMGViaApproxOrdersSetAO(DM dm, AO ao)
Set DM ordering.
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
MoFEMErrorCode DMMGViaApproxOrdersReplaceCoarseningIS(DM dm, IS *is_vec, int nb_elems, Mat A, int verb)
Replace coarsening IS in DM via approximation orders.
header of multi-grid solver for p- adaptivity
Post-process fields on refined mesh.
Implementation of solid shell prism element.
Implementation of solid shell prism element.
#define SSTR(x)
Convert number to string.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRQ_MOAB(a)
check error code of MoAB function
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMGetProblemFiniteElementLayout(DM dm, std::string fe_name, PetscLayout *layout)
Get finite elements layout in the problem.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
multi_index_container< boost::shared_ptr< DofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< DofEntity, UId, &DofEntity::getLocalUniqueId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< DofEntity, EntityHandle, &DofEntity::getEnt > > > > DofEntity_multiIndex
MultiIndex container keeps DofEntity.
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
implementation of Data Operators for Forces and Sources
PetscErrorCode pc_reset_mg(PC pc)
PetscErrorCode pc_setup_mg(PC pc)
PetscErrorCode pc_mg_set_last_level(PC pc, PetscInt levels, MPI_Comm *comms)
PetscErrorCode only_last_pc_reset_mg(PC pc)
static PetscErrorCode ierr
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
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
keeps basic data about problem
int nbLevels
number of multi-grid levels
std::map< EntityHandle, std::vector< EntityHandle > > elementsMap
MoFEMErrorCode generateReferenceElementMesh()
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
Set order of finite element.
PetscErrorCode operator()()
PetscErrorCode postProcess()
PetscErrorCode preProcess()
PetscErrorCode postProcess()
PetscErrorCode operator()()
PetscErrorCode preProcess()
PetscErrorCode getOptions()
get options from line command
PetscErrorCode createIsAtLevel(int kk, IS *is)
Set IS for levels.
PetscErrorCode buildProlongationOperator(PC pc, int verb=0)
PetscErrorCode destroyIsAtLevel(int kk, IS *is)
Destroy IS if internally created.
PetscErrorCode sortErrorVector(const double procent=0.33)
vector< int > vecOrderDofs
PetscErrorCode postProcFatPrims(const string &name)
Save post-processed data.
PetscErrorCode setOrder(int set, int up, double error, int max_order, int verb=0)
Set approximation order to entities.
PetscErrorCode printDisplacements()
Print displacements at selected points.
PetscErrorCode runKspAtLevel(DM dm, Mat Aij, Vec F, Vec D, bool add_level)
PetscErrorCode createIS(IS *is)
vector< Vec > solutionAtLevel
SolidShellPrismElement::CommonData & commonData
MoFEM::Interface & mField
PostProcFatPrismOnRefinedMesh postProcFatPrism
double procentErrorMax
max value at given precent
PetscErrorCode setUpOperators()
Set up operators for post-processing shell elements.
PetscErrorCode runKsp(DM dm, Mat Aij, Vec F, Vec D)
PetscErrorCode solveProblem(DM dm, Mat Aij, Vec F, Vec D, int nb_ref_cycles, int max_order, double procent=0.33, int verb=0)
Solve linear problem with p-adaptivity.
Calculate displacements at integration points.
Calculate rotation matrices to local (user) element coordinate system.
Get strain form displacements at integration points.
Get strain form local through thickness displacements at integration points.
Evaluate stress at integration points.
Calculate inverse of Jacobian.
Post process displacements on post-processing mesh.
Post process geometry (mesh nodal positions) on post-processing mesh.
Post process strain and stresses on post-processing mesh.