v0.14.0
Loading...
Searching...
No Matches
RunAdaptivity.cpp
Go to the documentation of this file.
1/** \file RunAdaptivity.cpp
2 * \ingroup solid_shell_prism_element
3 * \brief Implementation of solid shell prism element
4 *
5 */
6
7/*
8 * This file is part of MoFEM.
9 * MoFEM is free software: you can redistribute it and/or modify it under
10 * the terms of the GNU Lesser General Public License as published by the
11 * Free Software Foundation, either version 3 of the License, or (at your
12 * option) any later version.
13 *
14 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 * License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21
22#include <MoFEM.hpp>
24#include <boost/numeric/ublas/symmetric.hpp>
25
26using namespace MoFEM;
27#include <PostProcOnRefMesh.hpp>
28
29#ifdef WITH_ADOL_C
30// #define ADOLC_TAPELESS
31// #include <adolc/adtl.h>
32#include <adolc/adolc.h>
33#endif
34
36#include <DirichletBC.hpp>
38
39#include <ArcLengthTools.hpp>
41#include <RunAdaptivity.hpp>
42
43#if PETSC_VERSION_GE(3, 6, 0)
44#include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
45#include <petsc/private/pcmgimpl.h>
46#else
47#include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/
48#include <petsc-private/pcmgimpl.h>
49#endif
50
51namespace SolidShellModule {
52
53PetscErrorCode pc_setup_mg(PC pc);
54PetscErrorCode only_last_pc_reset_mg(PC pc);
55PetscErrorCode pc_mg_set_last_level(PC pc, PetscInt levels, MPI_Comm *comms);
56PetscErrorCode pc_reset_mg(PC pc);
57
58PetscErrorCode
60 IS *is) {
61 PetscFunctionBegin;
62 *is = mgLevels[kk];
63 PetscFunctionReturn(0);
64}
65
66PetscErrorCode
68 IS *is) {
69 PetscFunctionBegin;
70 PetscFunctionReturn(0);
71}
72
74 PetscFunctionBegin;
75 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
76 "MOFEM Shell Multi-Grid (Orders) pre-conditioner",
77 "none");
78
79 CHKERR PetscOptionsInt("-mofem_mg_verbose", "nb levels of multi-grid solver",
80 "", 0, &verboseLevel, PETSC_NULL);
81
82 ierr = PetscOptionsEnd();
83 CHKERRQ(ierr);
84
85 PetscFunctionReturn(0);
86}
87
88PetscErrorCode
90 PC pc, int verb) {
91 PetscFunctionBegin;
92 verb = verb > verboseLevel ? verb : verboseLevel;
93
94 MPI_Comm comm;
95 CHKERR PetscObjectGetComm((PetscObject)dM, &comm);
96
97 if (verb > 0) {
98 PetscPrintf(comm, "set MG levels %u\n", nbLevels);
99 }
100
101 vector<IS> is_vec(nbLevels + 1);
102 vector<int> is_glob_size(nbLevels + 1), is_loc_size(nbLevels + 1);
103
104 for (int kk = 0; kk < nbLevels; kk++) {
105
106 // get indices up to up to give approximation order
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]);
110
111 if (verb > 0) {
112 PetscSynchronizedPrintf(comm,
113 "Nb. dofs at level [ %d ] global %u local %d\n",
114 kk, is_glob_size[kk], is_loc_size[kk]);
115 }
116
117 // if no dofs on level kk finish here
118 if (is_glob_size[kk] == 0) {
119 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "no dofs at level");
120 }
121 }
122
123 // #if PETSC_VERSION_GE(3, 8, 0)
124 // #error "This is not working downgrade petsc-3.7.* or optimally petsc-3.6.*"
125 // #else
126 // CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
127 // CHKERRQ(ierr);
128 // #endif
129
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.*"
133#else
134 CHKERR PCMGSetGalerkin(pc, PETSC_TRUE);
135#endif
136
137 // CHKERR PCMGSetLevels(pc,nbLevels,NULL); CHKERRQ(ierr);
138 CHKERR pc_mg_set_last_level(pc, nbLevels, NULL);
139
140 CHKERR DMMGViaApproxOrdersReplaceCoarseningIS(dM, &is_vec[0], nbLevels, A, 1);
141
142 for (unsigned int kk = 0; kk < is_vec.size(); kk++) {
143 CHKERR destroyIsAtLevel(kk, &is_vec[kk]);
144 }
145
146 if (verb > 0) {
147 PetscSynchronizedFlush(comm, PETSC_STDOUT);
148 }
149
150 PetscFunctionReturn(0);
151}
152
154 PetscFunctionBegin;
155 PetscFunctionReturn(0);
156}
157
159 PetscFunctionBegin;
160 PetscFunctionReturn(0);
161}
162
164 PetscFunctionBegin;
165
166 EntityHandle ent = 0;
167 int order = 0;
168 int max_order = 0;
169
170 for (_IT_FENUMEREDDOF_BY_NAME_ROW_FOR_LOOP_(numeredEntFiniteElementPtr,
171 "DISPLACEMENT", dof)) {
172 if (dof->get()->getEntType() == MBVERTEX)
173 continue;
174 if (ent != dof->get()->getEnt()) {
175 ent = dof->get()->getEnt();
176 Range adj_prisms;
177 rval = mField.get_moab().get_adjacencies(&ent, 1, 3, false, adj_prisms);
179 vector<int> adj_orders(adj_prisms.size());
180 rval =
181 mField.get_moab().tag_get_data(thOrder, adj_prisms, &adj_orders[0]);
183 vector<int>::iterator max_it =
184 max_element(adj_orders.begin(), adj_orders.end());
185 max_order = *max_it;
186 if (order > max_order) {
187 SETERRQ2(
188 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
189 "order of element can not be bigger than order of entity %d > %d",
190 order, max_order);
191 }
192 rval = mField.get_moab().tag_set_data(thOrder, &ent, 1, &order);
194 if (order < max_order) {
195 rval = mField.get_moab().tag_set_data(thOrder, &ent, 1, &max_order);
197 }
198 }
199 }
200
201 PetscFunctionReturn(0);
202}
203
205 PetscFunctionBegin;
206 {
207 double def_val = 0;
208 rval = mField.get_moab().tag_get_handle(
209 "ENERGY_OR_ENERGY_ERROR", 1, MB_TYPE_DOUBLE, thError,
210 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
212 }
213 PetscFunctionReturn(0);
214}
215
217 PetscFunctionBegin;
218 PetscFunctionReturn(0);
219}
220
222 PetscFunctionBegin;
223 EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
224 int order;
225 rval = mField.get_moab().tag_get_data(thOrder, &fe_ent, 1, &order);
227
228 if (order == maxOrder)
229 PetscFunctionReturn(0);
230
231 if (eRror > 0) {
232 double elem_error;
233 rval = mField.get_moab().tag_get_data(thError, &fe_ent, 1, &elem_error);
235 if (elem_error < eRror && order > 2) {
236 order -= 1;
237 rval = mField.get_moab().tag_set_data(thOrder, &fe_ent, 1, &order);
239 }
240 }
241
242 if (sEt > 0) {
243 order = sEt;
244 } else {
245 if (order < maxOrder) {
246 order += uP;
247 }
248 }
249
250 rval = mField.get_moab().tag_set_data(thOrder, &fe_ent, 1, &order);
252
253 PetscFunctionReturn(0);
254}
255
256PetscErrorCode RunAdaptivity::setOrder(int set, int up, double error,
257 int max_order, int verb) {
258 PetscFunctionBegin;
259 int def_val = 1;
260 rval = mField.get_moab().tag_get_handle("ADAPT_ORDER", 1, MB_TYPE_INTEGER,
261 thOrder, MB_TAG_CREAT | MB_TAG_SPARSE,
262 &def_val);
264 SetOrderToElement fe_order_elements(mField, thOrder, set, up, error,
265 max_order);
266 CHKERR DMoFEMLoopFiniteElements(dM, "SHELL", &fe_order_elements);
268 Range prisms;
269
270 rval = mField.get_moab().get_entities_by_handle(meshset, prisms);
272 Tag th;
273 rval = mField.get_moab().tag_get_handle("ADAPT_ORDER", th);
275 ParallelComm *pcomm =
276 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
277 // rval = pcomm->reduce_tags(th,MPI_SUM,prisms);
278 rval = pcomm->exchange_tags(th, prisms);
279 if (verb > 1) {
280 vector<Tag> tags;
281 tags.push_back(th);
282 tags.push_back(pcomm->part_tag());
283 Tag th_global_id;
284 rval = mField.get_moab().tag_get_handle(GLOBAL_ID_TAG_NAME, th_global_id);
286 tags.push_back(th_global_id);
287 for (int rr = 0; rr < mField.get_comm_size(); rr++) {
288 if (mField.get_comm_rank() == rr) {
289 ostringstream o1;
290 o1 << "test_" << rr << ".vtk";
291 mField.get_moab().write_file(o1.str().c_str(), "VTK", "", &meshset, 1,
292 &tags[0], tags.size());
293 }
294 }
295 }
296 // CHKERR PetscBarrier(PETSC_NULL); CHKERRQ(ierr);
297 SetOrderToEntities fe_order_entities(mField, thOrder);
298 CHKERR DMoFEMLoopFiniteElements(dM, "SHELL", &fe_order_entities);
299 PetscFunctionReturn(0);
300}
301
302PetscErrorCode RunAdaptivity::createIS(IS *is) {
303 PetscFunctionBegin;
304 const MoFEM::Problem *problem_ptr;
305 CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
306 EntityHandle ent = 0;
307 int ent_order;
308 vecOrderDofs.clear();
309 for (_IT_NUMEREDDOF_ROW_BY_OWNPROC_FOR_LOOP_(problem_ptr,
310 mField.get_comm_rank(), dof)) {
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());
315 } else {
316 if (ent != dof->get()->getEnt()) {
317 ent = dof->get()->getEnt();
318 rval = mField.get_moab().tag_get_data(thOrder, &ent, 1, &ent_order);
320 }
321 if (dof->get()->getDofOrder() <= ent_order) {
322 vecOrderDofs.push_back(dof->get()->getPetscGlobalDofIdx());
323 }
324 }
325 }
326 CHKERR ISCreateGeneral(PETSC_COMM_WORLD, vecOrderDofs.size(),
327 &*vecOrderDofs.begin(), PETSC_COPY_VALUES, is);
328 PetscFunctionReturn(0);
329}
330
331PetscErrorCode RunAdaptivity::runKsp(DM dm, Mat Aij, Vec F, Vec D) {
332 PetscFunctionBegin;
333 PetscPrintf(PETSC_COMM_WORLD, "SetUp sOlver ... ");
334 CHKERR KSPSetDM(sOlver, dm);
335 CHKERR KSPSetOperators(sOlver, Aij, Aij);
336 CHKERR KSPSetDMActive(sOlver, PETSC_FALSE);
337 if (1) {
338 // from PETSc example ex42.c
339 PetscBool same = PETSC_FALSE;
340 CHKERR KSPGetPC(sOlver, &pC);
341 PetscErrorCode (*org_pc_setup_mg)(PC pc) = pC->ops->setup;
342 pC->ops->setup = pc_setup_mg;
343 PetscObjectTypeCompare((PetscObject)pC, PCMG, &same);
344 if (same) {
345 PetscPrintf(PETSC_COMM_WORLD, "multigrid ..\n");
346 int nb_levels = mgLevels.size();
347 AO ao;
348 CHKERR AOCreateMappingIS(mgLevels.back(), PETSC_NULL, &ao);
351 pc_ctx.nbLevels = nb_levels;
353 CHKERR AODestroy(&ao);
354 } else {
355 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
356 "Expected MG pre-conditioner");
357 }
358 CHKERR KSPSetInitialGuessKnoll(sOlver, PETSC_FALSE);
359 CHKERR KSPSetInitialGuessNonzero(sOlver, PETSC_TRUE);
360 PetscPrintf(PETSC_COMM_WORLD, "Solve problem ..\n");
361 CHKERR KSPSetUp(sOlver);
362 CHKERR KSPSolve(sOlver, F, D);
363 PetscErrorCode (*org_reset)(PC pc) = pC->ops->reset;
364 pC->ops->reset = only_last_pc_reset_mg;
365 CHKERR KSPReset(sOlver);
366 pC->ops->reset = org_reset;
367 pC->ops->setup = org_pc_setup_mg;
368 }
369 PetscFunctionReturn(0);
370}
371
372PetscErrorCode RunAdaptivity::runKspAtLevel(DM dm, Mat Aij, Vec F, Vec D,
373 bool add_level) {
374 PetscFunctionBegin;
375 PetscPrintf(PETSC_COMM_WORLD, "Get level sub matrix ...\n");
376 Mat sub_Aij;
377 CHKERR MatGetSubMatrix(Aij, mgLevels.back(), mgLevels.back(),
378 MAT_INITIAL_MATRIX, &sub_Aij);
379 // {
380 // CHKERR MatView(sub_Aij,PETSC_VIEWER_DRAW_WORLD); CHKERRQ(ierr);
381 // std::string wait;
382 // std::cin >> wait;
383 // }
384 Vec sub_F;
385 CHKERR VecGetSubVector(F, mgLevels.back(), &sub_F);
386 Vec sub_D;
387 CHKERR VecGetSubVector(D, mgLevels.back(), &sub_D);
388 CHKERR runKsp(dm, sub_Aij, sub_F, sub_D);
389 CHKERR VecRestoreSubVector(F, mgLevels.back(), &sub_F);
390 CHKERR VecRestoreSubVector(D, mgLevels.back(), &sub_D);
391 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
392 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
393 CHKERR MatDestroy(&sub_Aij);
394 PetscFunctionReturn(0);
395}
396
397PetscErrorCode RunAdaptivity::sortErrorVector(const double procent) {
398 PetscFunctionBegin;
399 int size;
400 CHKERR VecGetSize(commonData.vecErrorOnElements, &size);
401 int loc_size;
402 CHKERR VecGetLocalSize(commonData.vecErrorOnElements, &loc_size);
403 double *a0;
404 CHKERR VecGetArray(commonData.vecErrorOnElements, &a0);
405 for (int ii = 0; ii < loc_size; ii++) {
406 a0[ii] = 1 / a0[ii];
407 }
408 CHKERR PetscSortReal(loc_size, a0);
409 int procent_size = ceil(size * procent);
410 procent_size = procent_size > loc_size ? loc_size : procent_size;
411 Vec short_sorted;
412 CHKERR VecCreateMPIWithArray(mField.get_comm(), 1, procent_size,
413 PETSC_DETERMINE, a0, &short_sorted);
414 if (mField.get_comm_size() > 1) {
415 VecScatter ctx;
416 Vec all_in_zero;
417 CHKERR VecScatterCreateToZero(short_sorted, &ctx, &all_in_zero);
418 CHKERR VecScatterBegin(ctx, short_sorted, all_in_zero, INSERT_VALUES,
419 SCATTER_FORWARD);
420 CHKERR VecScatterEnd(ctx, short_sorted, all_in_zero, INSERT_VALUES,
421 SCATTER_FORWARD);
422 if (mField.get_comm_rank() == 0) {
423 double *a1;
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;
430 procentErrorMax = a1[pos];
432 CHKERR VecRestoreArray(all_in_zero, &a1);
433 }
434 CHKERR VecDestroy(&all_in_zero);
435 CHKERR VecScatterDestroy(&ctx);
436 CHKERR VecDestroy(&short_sorted);
437 Vec exchange_value;
438 int one[] = {0};
439 bool rank0 = (unsigned int)mField.get_comm_rank() == 0;
440 CHKERR VecCreateGhostWithArray(mField.get_comm(), rank0 ? 1 : 0, 1,
441 rank0 ? 0 : 1, one, &procentErrorMax,
442 &exchange_value);
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++) {
447 a0[ii] = 1 / a0[ii];
448 }
449 } else {
450 procentErrorMax = a0[(unsigned int)ceil(size * procent)];
452 }
453 CHKERR VecDestroy(&short_sorted);
454 CHKERR VecRestoreArray(commonData.vecErrorOnElements, &a0);
455 PetscFunctionReturn(0);
456}
457
458PetscErrorCode RunAdaptivity::solveProblem(DM dm, Mat Aij, Vec F, Vec D,
459 int nb_ref_cycles, int max_order,
460 double procent, int verb) {
461 PetscFunctionBegin;
462
463 const int max_levels = 40;
464 isLevels.resize(max_levels, PETSC_NULL);
465 solutionAtLevel.resize(max_levels, PETSC_NULL);
466
467 CHKERR KSPCreate(PETSC_COMM_WORLD, &sOlver);
468 CHKERR KSPSetFromOptions(sOlver);
469
470 PetscLayout layout;
472 int size;
473 CHKERR PetscLayoutGetLocalSize(layout, &size);
474 CHKERR VecCreateMPI(PETSC_COMM_WORLD, size, PETSC_DETERMINE,
476 CHKERR PetscLayoutDestroy(&layout);
477
478 // Solve level 0
479 CHKERR setOrder(1, 0, 0, max_order);
481 mgLevels.push_back(isLevels[0]);
482
483 // Solve level 1
484 CHKERR setOrder(0, 1, 0, max_order);
486 mgLevels.push_back(isLevels[1]);
487 CHKERR runKspAtLevel(dm, Aij, F, D, true);
488 CHKERR VecDuplicate(D, &solutionAtLevel[1]);
489 CHKERR VecCopy(D, solutionAtLevel[1]);
490
491 commonData.saveOnTag = false;
492 CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_REVERSE);
494 double energy;
495 CHKERR VecSum(commonData.vecErrorOnElements, &energy);
496 PetscPrintf(PETSC_COMM_WORLD, "Level %d Energy = %12.9e\n", 0, energy);
497
498 if (verb > 0) {
500 CHKERR postProcFatPrims("solid_" + SSTR(0) + ".h5m");
501 }
502
503 Vec E;
504 CHKERR VecDuplicate(D, &E);
505
506 int rr = 1;
507 for (;; rr++) {
508
509 double max_error = 0;
510 if (rr % 2 == 0) {
511 double vnorm;
512 int is_size;
513 {
514 CHKERR VecWAXPY(E, -1, solutionAtLevel[rr - 1], solutionAtLevel[rr]);
515 CHKERR VecGhostUpdateBegin(E, INSERT_VALUES, SCATTER_FORWARD);
516 CHKERR VecGhostUpdateEnd(E, INSERT_VALUES, SCATTER_FORWARD);
517 CHKERR VecNorm(E, NORM_2, &vnorm);
518 CHKERR DMoFEMMeshToLocalVector(dM, E, INSERT_VALUES, SCATTER_REVERSE);
519 commonData.saveOnTag = true;
521 CHKERR ISGetSize(isLevels[rr], &is_size);
522 CHKERR VecMax(commonData.vecErrorOnElements, PETSC_NULL, &max_error);
523 CHKERR sortErrorVector(procent);
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",
527 rr, is_size, vnorm, procentErrorMax, max_error);
528 }
529 {
530 CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_REVERSE);
531 double vnorm;
532 CHKERR VecNorm(D, NORM_2, &vnorm);
533 commonData.saveOnTag = false;
535 double energy = 0;
536 CHKERR VecSum(commonData.vecErrorOnElements, &energy);
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);
540 }
541 if (verb > 0) {
543 CHKERR postProcFatPrims("solid_" + SSTR(rr + 1) + ".h5m");
544 }
545 }
546 if (rr == nb_ref_cycles) {
547 break;
548 }
549
550 if (rr % 2 == 0) {
551 CHKERR setOrder(0, 0, procentErrorMax, max_order);
552 } else {
553 CHKERR setOrder(0, 1, 0, max_order);
554 }
555
556 CHKERR createIS(&isLevels[rr + 1]);
557 mgLevels.push_back(isLevels[rr + 1]);
558 if (mgLevels.size() > 1) {
559 int s1, s2;
560 CHKERR ISGetSize(mgLevels.back(), &s1);
561 CHKERR ISGetSize(mgLevels[mgLevels.size() - 2], &s2);
562 PetscPrintf(PETSC_COMM_WORLD,
563 "Last level size = %d, previous level size = %d\n", s1, s2);
564 if (s1 == s2) {
565 IS is = mgLevels.back();
566 CHKERR ISDestroy(&mgLevels[mgLevels.size() - 2]);
567 mgLevels.pop_back();
568 mgLevels.back() = is;
569 isLevels[rr] = PETSC_NULL;
570 }
571 }
572 {
573 Vec sub_d, sub_b;
574 VecGetSubVector(D, mgLevels[mgLevels.size() - 2], &sub_d);
575 VecGetSubVector(solutionAtLevel[rr], mgLevels[mgLevels.size() - 2],
576 &sub_b);
577 CHKERR VecCopy(sub_b, sub_d);
578 VecRestoreSubVector(D, mgLevels[mgLevels.size() - 2], &sub_d);
579 VecRestoreSubVector(solutionAtLevel[rr], mgLevels[mgLevels.size() - 2],
580 &sub_b);
581 }
582 CHKERR runKspAtLevel(dm, Aij, F, D, rr % 2 == 0);
583 CHKERR VecDuplicate(D, &solutionAtLevel[rr + 1]);
584 {
585 Vec sub_d, sub_b;
586 VecGetSubVector(D, mgLevels.back(), &sub_d);
587 VecGetSubVector(solutionAtLevel[rr + 1], mgLevels.back(), &sub_b);
588 CHKERR VecCopy(sub_d, sub_b);
589 VecRestoreSubVector(D, mgLevels.back(), &sub_d);
590 VecRestoreSubVector(solutionAtLevel[rr + 1], mgLevels.back(), &sub_b);
591 }
592
593 if (rr % 2 != 0) {
594 mgLevels.pop_back();
595 }
596 }
597
598 for (int ll = 0; ll < max_levels; ll++) {
599 // PetscPrintf(PETSC_COMM_WORLD,"Destroy %d of %d\n",ll,rr);
600 if (isLevels[ll]) {
601 // cerr << "destroy " << ll << endl;
602 CHKERR ISDestroy(&isLevels[ll]);
603 }
604 if (solutionAtLevel[ll]) {
605 CHKERR VecDestroy(&solutionAtLevel[ll]);
606 }
607 }
608
610 CHKERR VecDestroy(&E);
611
612 CHKERR KSPDestroy(&sOlver);
613
614 PetscFunctionReturn(0);
615}
616
618 PetscFunctionBegin;
647 commonData));
651 commonData));
653 PetscFunctionReturn(0);
654}
655
656PetscErrorCode RunAdaptivity::postProcFatPrims(const string &name) {
657 PetscFunctionBegin;
659 // Save data on mesh
660 rval = postProcFatPrism.postProcMesh.write_file(name.c_str(), "MOAB",
661 "PARALLEL=WRITE_PART");
663 PetscFunctionReturn(0);
664}
665
667 PetscFunctionBegin;
668 const DofEntity_multiIndex *dofs_ptr;
669 CHKERR mField.get_dofs(&dofs_ptr);
670 for (_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(mField, "RESULTS", it)) {
671 EntityHandle meshset = it->getMeshset();
672 Range nodes;
673 rval =
674 mField.get_moab().get_entities_by_type(meshset, MBVERTEX, nodes, true);
676 Range edges;
677 rval = mField.get_moab().get_adjacencies(nodes, 1, false, edges,
678 moab::Interface::UNION);
680 Range nodes_faces;
681 rval = mField.get_moab().get_adjacencies(nodes, 2, false, nodes_faces,
682 moab::Interface::UNION);
684 Range tris = nodes_faces.subset_by_type(MBTRI);
685 Range tris_nodes;
686 rval = mField.get_moab().get_connectivity(tris, tris_nodes, true);
688 Range faces;
689 rval = mField.get_moab().get_adjacencies(edges, 2, false, faces,
690 moab::Interface::UNION);
692 Range quads;
693 quads = faces.subset_by_type(MBQUAD);
694 Range edges_nodes;
695 rval = mField.get_moab().get_connectivity(edges, edges_nodes, true);
697 Range quad_nodes;
698 rval = mField.get_moab().get_connectivity(quads, quad_nodes, true);
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;
704 dit = dofs_ptr->get<Composite_Name_And_Ent_mi_tag>().lower_bound(
705 boost::make_tuple("DISPLACEMENT", *nit));
706 hi_dit = dofs_ptr->get<Composite_Name_And_Ent_mi_tag>().upper_bound(
707 boost::make_tuple("DISPLACEMENT", *nit));
708 for (; dit != hi_dit; dit++) {
709 EntityHandle ent = dit->get()->getEnt();
710 double coords[3];
711 rval = mField.get_moab().get_coords(&*nit, 1, coords);
713 PetscSynchronizedPrintf(
714 PETSC_COMM_WORLD,
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],
718 coords[1], coords[2], ent, mField.get_comm_rank());
719 }
720 }
721 }
722 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
723 PetscFunctionReturn(0);
724}
725
726} // namespace SolidShellModule
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
Definition definitions.h:31
#define CHKERR
Inline error check.
constexpr int order
@ F
PetscErrorCode DMMoFEMGetProblemFiniteElementLayout(DM dm, std::string fe_name, PetscLayout *layout)
Get finite elements layout in the problem.
Definition DMMoFEM.cpp:459
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:509
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:412
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
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.
constexpr double a0
constexpr double a1
double D
constexpr int nb_levels
Definition level_set.cpp:58
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode pc_reset_mg(PC pc)
PetscErrorCode pc_setup_mg(PC pc)
Definition HackMG.cpp:227
PetscErrorCode pc_mg_set_last_level(PC pc, PetscInt levels, MPI_Comm *comms)
Definition HackMG.cpp:73
PetscErrorCode only_last_pc_reset_mg(PC pc)
Definition HackMG.cpp:21
static PetscErrorCode ierr
Definition HackMG.cpp:17
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
PetscErrorCode getOptions()
get options from line command
PetscErrorCode createIsAtLevel(int kk, IS *is)
Set IS for levels.
PetscErrorCode destroyIsAtLevel(int kk, IS *is)
Destroy IS if internally created.
PetscErrorCode sortErrorVector(const double procent=0.33)
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)
SolidShellPrismElement::CommonData & commonData
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 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.
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.