v0.15.0
Loading...
Searching...
No Matches
Smoother.hpp
Go to the documentation of this file.
1
2
3#ifndef __SMOOTHER_HPP__
4#define __SMOOTHER_HPP__
5
6#ifndef WITH_ADOL_C
7#error "MoFEM need to be compiled with ADOL-C"
8#endif
9
10struct Smoother {
11
12 struct SmootherBlockData {
13
15 Vec frontF;
18
20 : sTabilised(false), frontF(PETSC_NULLPTR), tangentFrontF(PETSC_NULLPTR),
21 ownVectors(false) {
22 ierr = getOptions();
23 CHKERRABORT(PETSC_COMM_SELF, ierr);
24 }
25
26 MoFEMErrorCode getOptions() {
28 PetscOptionsBegin(PETSC_COMM_WORLD, "",
29 "Get stabilisation element options", "none");
31 PetscBool smoothing_on = sTabilised ? PETSC_TRUE : PETSC_FALSE;
32 CHKERR PetscOptionsBool("-smoothing_stabilise",
33 "all nodes controlled by smoothing element", "",
34 smoothing_on, &smoothing_on, PETSC_NULLPTR);
35 sTabilised = (smoothing_on == PETSC_TRUE) ? true : false;
36 PetscOptionsEnd();
38 }
39
41 if (ownVectors) {
42 ierr = VecDestroy(&frontF);
43 CHKERRABORT(PETSC_COMM_WORLD, ierr);
44 ierr = VecDestroy(&tangentFrontF);
45 CHKERRABORT(PETSC_COMM_WORLD, ierr);
46 }
47 }
48 };
50
51 std::map<int, NonlinearElasticElement::BlockData> setOfBlocks;
53
55
57
60 smootherData(smoother_data) {}
61
62 MoFEMErrorCode preProcess() {
64
65 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
66
67 if (A != PETSC_NULLPTR)
68 snes_B = A;
69
70 if (F != PETSC_NULLPTR)
71 snes_f = F;
72 switch (snes_ctx) {
74 if (smootherData.frontF) {
75 CHKERR VecZeroEntries(smootherData.frontF);
76 CHKERR VecGhostUpdateBegin(smootherData.frontF, INSERT_VALUES,
77 SCATTER_FORWARD);
78 CHKERR VecGhostUpdateEnd(smootherData.frontF, INSERT_VALUES,
79 SCATTER_FORWARD);
80 }
81 } break;
82 default:
83 break;
84 }
85
87 }
88
89 MoFEMErrorCode postProcess() {
91
92 switch (snes_ctx) {
94 if (smootherData.frontF) {
95 CHKERR VecAssemblyBegin(smootherData.frontF);
96 CHKERR VecAssemblyEnd(smootherData.frontF);
97 CHKERR VecGhostUpdateBegin(smootherData.frontF, ADD_VALUES,
98 SCATTER_REVERSE);
99 CHKERR VecGhostUpdateEnd(smootherData.frontF, ADD_VALUES,
100 SCATTER_REVERSE);
101 CHKERR VecGhostUpdateBegin(smootherData.frontF, INSERT_VALUES,
102 SCATTER_FORWARD);
103 CHKERR VecGhostUpdateEnd(smootherData.frontF, INSERT_VALUES,
104 SCATTER_FORWARD);
105 }
106 break;
107 default:
108 break;
109 }
110 }
111
112 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
114 }
115 };
116
117 boost::shared_ptr<MyVolumeFE> feRhsPtr;
118 boost::shared_ptr<MyVolumeFE> feLhsPtr;
119
120 MyVolumeFE &feRhs; ///< calculate right hand side for tetrahedral elements
121 MyVolumeFE &getLoopFeRhs() { return feRhs; } ///< get rhs volume element
122 MyVolumeFE &feLhs; //< calculate left hand side for tetrahedral elements
123 MyVolumeFE &getLoopFeLhs() { return feLhs; } ///< get lhs volume element
124
126 : feRhsPtr(new MyVolumeFE(m_field, smootherData)),
128 feLhs(*feLhsPtr) {}
129
130 struct OpJacobianSmoother
132
136 int tag, bool jacobian)
138 field_name, data, common_data, tag, jacobian, false, false) {}
139
140 MoFEMErrorCode calculateStress(const int gg) {
142
143 CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
145
146 commonData.sTress[gg].resize(3, 3, false);
147 for (int dd1 = 0; dd1 < 3; dd1++) {
148 for (int dd2 = 0; dd2 < 3; dd2++) {
149 dAta.materialAdoublePtr->P(dd1, dd2) >>=
150 (commonData.sTress[gg])(dd1, dd2);
151 }
152 }
153
155 }
156 };
157
158 struct OpRhsSmoother : public NonlinearElasticElement::OpRhsPiolaKirchhoff {
159
161
162 OpRhsSmoother(const std::string field_name,
165 SmootherBlockData &smoother_data)
167 common_data),
168 smootherData(smoother_data) {}
169
170 ublas::vector<int> frontIndices;
171
172 MoFEMErrorCode aSemble(int row_side, EntityType row_type,
173 EntitiesFieldData::EntData &row_data) {
175
176 int nb_dofs = row_data.getIndices().size();
177 int *indices_ptr = &row_data.getIndices()[0];
178
179 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
180 iNdices.resize(nb_dofs, false);
181 noalias(iNdices) = row_data.getIndices();
183 indices_ptr = &iNdices[0];
184 }
185 frontIndices.resize(nb_dofs, false);
186 noalias(frontIndices) = row_data.getIndices();
187 VectorDofs &dofs = row_data.getFieldDofs();
188 VectorDofs::iterator dit = dofs.begin();
189 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
190 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) !=
192 iNdices[ii] = -1;
193 } else {
194 frontIndices[ii] = -1;
195 }
196 }
197 if (smootherData.frontF) {
198 CHKERR VecSetValues(smootherData.frontF, nb_dofs, &frontIndices[0],
199 &nf[0], ADD_VALUES);
200 }
201 }
202
203 CHKERR VecSetOption(getFEMethod()->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
204 PETSC_TRUE);
205 CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
206 ADD_VALUES);
207
209 }
210 };
211
212 struct OpLhsSmoother
214
217
219 const std::string vel_field, const std::string field_name,
222 SmootherBlockData &smoother_data,
223 const std::string
224 crack_area_tangent_constrains // = "LAMBDA_CRACK_TANGENT_CONSTRAIN"
225 )
227 data, common_data),
228 smootherData(smoother_data),
229 fieldCrackAreaTangentConstrains(crack_area_tangent_constrains) {}
230
231 ublas::vector<int> rowFrontIndices;
232
233 MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
234 EntityType col_type,
235 EntitiesFieldData::EntData &row_data,
236 EntitiesFieldData::EntData &col_data) {
238
239 int nb_row = row_data.getIndices().size();
240 int nb_col = col_data.getIndices().size();
241 int *row_indices_ptr = &row_data.getIndices()[0];
242 int *col_indices_ptr = &col_data.getIndices()[0];
243
244 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
245 rowIndices.resize(nb_row, false);
246 noalias(rowIndices) = row_data.getIndices();
248 row_indices_ptr = &rowIndices[0];
249 }
250 rowFrontIndices.resize(nb_row, false);
251 noalias(rowFrontIndices) = row_data.getIndices();
252 VectorDofs &dofs = row_data.getFieldDofs();
253 VectorDofs::iterator dit = dofs.begin();
254 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
255 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) !=
257 rowIndices[ii] = -1;
258 } else {
259 rowFrontIndices[ii] = -1;
260 }
261 }
262 }
263
264 CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr,
265 nb_col, col_indices_ptr, &k(0, 0), ADD_VALUES);
266
268
269 const auto bit_number_for_crack_area_tangent_constrain =
271 const auto bit_number_for_mesh_position =
272 getFEMethod()->getFieldBitNumber("MESH_NODE_POSITIONS");
273
274 // get tangent vector array
275 double *f_tangent_front_mesh_array;
277 &f_tangent_front_mesh_array);
278
279 auto row_dofs = getFEMethod()->getRowDofsPtr();
280
281 // iterate nodes on tet
282 for (int nn = 0; nn < 4; nn++) {
283
284 // get indices with Lagrange multiplier at node nn
285 auto dit = row_dofs->get<Unique_mi_tag>().lower_bound(
286 DofEntity::getLoFieldEntityUId(
287 bit_number_for_crack_area_tangent_constrain, getConn()[nn]));
288 auto hi_dit = row_dofs->get<Unique_mi_tag>().upper_bound(
289 DofEntity::getHiFieldEntityUId(
290 bit_number_for_crack_area_tangent_constrain, getConn()[nn]));
291
292 // continue if Lagrange are on element
293 if (std::distance(dit, hi_dit) > 0) {
294
295 // get mesh node positions at node nn
296 auto diit = row_dofs->get<Unique_mi_tag>().lower_bound(
297 DofEntity::getLoFieldEntityUId(bit_number_for_mesh_position,
298 getConn()[nn]));
299
300 auto hi_diit = row_dofs->get<Unique_mi_tag>().upper_bound(
301 DofEntity::getHiFieldEntityUId(bit_number_for_mesh_position,
302 getConn()[nn]));
303
304 // iterate over dofs on node nn
305 for (; diit != hi_diit; diit++) {
306 // iterate overt dofs in element column
307 for (int ddd = 0; ddd < nb_col; ddd++) {
308 // check consistency, node has to be at crack front
309 if (rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()] !=
310 diit->get()->getPetscGlobalDofIdx()) {
311 SETERRQ(
312 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
313 "data inconsistency %d != %d",
314 rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()],
315 diit->get()->getPetscGlobalDofIdx());
316 }
317 // dof is not on this partition
318 if (diit->get()->getPetscLocalDofIdx() == -1) {
319 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
320 "data inconsistency");
321 }
322 double g =
323 f_tangent_front_mesh_array[diit->get()
324 ->getPetscLocalDofIdx()] *
325 k(3 * nn + diit->get()->getDofCoeffIdx(), ddd);
326 int lambda_idx = dit->get()->getPetscGlobalDofIdx();
327 CHKERR MatSetValues(getFEMethod()->snes_B, 1, &lambda_idx, 1,
328 &col_indices_ptr[ddd], &g, ADD_VALUES);
329 }
330 }
331 }
332 }
333 CHKERR VecRestoreArray(smootherData.tangentFrontF,
334 &f_tangent_front_mesh_array);
335 }
336
338 }
339 };
340};
341
342#endif //__SMOOTHER_HPP__
static PetscErrorCode ierr
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr auto field_name
constexpr double g
unsigned int getFieldBitNumber(std::string field_name) const
Deprecated interface functions.
auto getRowDofsPtr() const
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Vec & snes_f
residual
Mat & snes_B
preconditioner of jacobian matrix
SNESContext snes_ctx
data for calculation heat conductivity and heat capacity elements
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr
common data used by volume elements
OpJacobianPiolaKirchhoffStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale, bool field_disp)
Construct operator to calculate Piola-Kirchhoff stress or its derivatives over gradient deformation.
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
structure grouping operators and data used for calculation of nonlinear elastic element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition Smoother.hpp:62
MyVolumeFE(MoFEM::Interface &m_field, SmootherBlockData &smoother_data)
Definition Smoother.hpp:58
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition Smoother.hpp:89
SmootherBlockData & smootherData
Definition Smoother.hpp:56
OpJacobianSmoother(const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, int tag, bool jacobian)
Definition Smoother.hpp:133
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
Definition Smoother.hpp:140
ublas::vector< int > rowFrontIndices
Definition Smoother.hpp:231
OpLhsSmoother(const std::string vel_field, const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, SmootherBlockData &smoother_data, const std::string crack_area_tangent_constrains)
Definition Smoother.hpp:218
const std::string fieldCrackAreaTangentConstrains
Definition Smoother.hpp:216
SmootherBlockData & smootherData
Definition Smoother.hpp:215
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition Smoother.hpp:233
ublas::vector< int > frontIndices
Definition Smoother.hpp:170
SmootherBlockData & smootherData
Definition Smoother.hpp:160
OpRhsSmoother(const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, SmootherBlockData &smoother_data)
Definition Smoother.hpp:162
MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition Smoother.hpp:172
MoFEMErrorCode getOptions()
Definition Smoother.hpp:26
std::map< int, NonlinearElasticElement::BlockData > setOfBlocks
Definition Smoother.hpp:51
MyVolumeFE & feLhs
Definition Smoother.hpp:122
SmootherBlockData smootherData
Definition Smoother.hpp:49
MyVolumeFE & getLoopFeLhs()
get lhs volume element
Definition Smoother.hpp:123
Smoother(MoFEM::Interface &m_field)
Definition Smoother.hpp:125
NonlinearElasticElement::CommonData commonData
Definition Smoother.hpp:52
boost::shared_ptr< MyVolumeFE > feLhsPtr
Definition Smoother.hpp:118
MyVolumeFE & getLoopFeRhs()
get rhs volume element
Definition Smoother.hpp:121
boost::shared_ptr< MyVolumeFE > feRhsPtr
Definition Smoother.hpp:117
MyVolumeFE & feRhs
calculate right hand side for tetrahedral elements
Definition Smoother.hpp:120