v0.14.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
13
15 Vec frontF;
18
20 : sTabilised(false), frontF(PETSC_NULL), tangentFrontF(PETSC_NULL),
21 ownVectors(false) {
22 ierr = getOptions();
23 CHKERRABORT(PETSC_COMM_SELF, ierr);
24 }
25
26 MoFEMErrorCode getOptions() {
28 ierr = 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_NULL);
35 sTabilised = (smoothing_on == PETSC_TRUE) ? true : false;
36 ierr = PetscOptionsEnd();
39 }
40
42 if (ownVectors) {
43 ierr = VecDestroy(&frontF);
44 CHKERRABORT(PETSC_COMM_WORLD, ierr);
45 ierr = VecDestroy(&tangentFrontF);
46 CHKERRABORT(PETSC_COMM_WORLD, ierr);
47 }
48 }
49 };
51
52 std::map<int, NonlinearElasticElement::BlockData> setOfBlocks;
54
56
58
61 smootherData(smoother_data) {}
62
63 MoFEMErrorCode preProcess() {
65
66 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
67
68 if (A != PETSC_NULL)
69 snes_B = A;
70
71 if (F != PETSC_NULL)
72 snes_f = F;
73 switch (snes_ctx) {
75 if (smootherData.frontF) {
76 CHKERR VecZeroEntries(smootherData.frontF);
77 CHKERR VecGhostUpdateBegin(smootherData.frontF, INSERT_VALUES,
78 SCATTER_FORWARD);
79 CHKERR VecGhostUpdateEnd(smootherData.frontF, INSERT_VALUES,
80 SCATTER_FORWARD);
81 }
82 } break;
83 default:
84 break;
85 }
86
88 }
89
90 MoFEMErrorCode postProcess() {
92
93 switch (snes_ctx) {
95 if (smootherData.frontF) {
96 CHKERR VecAssemblyBegin(smootherData.frontF);
97 CHKERR VecAssemblyEnd(smootherData.frontF);
98 CHKERR VecGhostUpdateBegin(smootherData.frontF, ADD_VALUES,
99 SCATTER_REVERSE);
100 CHKERR VecGhostUpdateEnd(smootherData.frontF, ADD_VALUES,
101 SCATTER_REVERSE);
102 CHKERR VecGhostUpdateBegin(smootherData.frontF, INSERT_VALUES,
103 SCATTER_FORWARD);
104 CHKERR VecGhostUpdateEnd(smootherData.frontF, INSERT_VALUES,
105 SCATTER_FORWARD);
106 }
107 break;
108 default:
109 break;
110 }
111 }
112
113 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
115 }
116 };
117
118 boost::shared_ptr<MyVolumeFE> feRhsPtr;
119 boost::shared_ptr<MyVolumeFE> feLhsPtr;
120
121 MyVolumeFE &feRhs; ///< calculate right hand side for tetrahedral elements
122 MyVolumeFE &getLoopFeRhs() { return feRhs; } ///< get rhs volume element
123 MyVolumeFE &feLhs; //< calculate left hand side for tetrahedral elements
124 MyVolumeFE &getLoopFeLhs() { return feLhs; } ///< get lhs volume element
125
127 : feRhsPtr(new MyVolumeFE(m_field, smootherData)),
129 feLhs(*feLhsPtr) {}
130
133
137 int tag, bool jacobian)
139 field_name, data, common_data, tag, jacobian, false, false) {}
140
141 MoFEMErrorCode calculateStress(const int gg) {
143
144 CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
146
147 commonData.sTress[gg].resize(3, 3, false);
148 for (int dd1 = 0; dd1 < 3; dd1++) {
149 for (int dd2 = 0; dd2 < 3; dd2++) {
150 dAta.materialAdoublePtr->P(dd1, dd2) >>=
151 (commonData.sTress[gg])(dd1, dd2);
152 }
153 }
154
156 }
157 };
158
160
162
163 OpRhsSmoother(const std::string field_name,
166 SmootherBlockData &smoother_data)
168 common_data),
169 smootherData(smoother_data) {}
170
171 ublas::vector<int> frontIndices;
172
173 MoFEMErrorCode aSemble(int row_side, EntityType row_type,
174 EntitiesFieldData::EntData &row_data) {
176
177 int nb_dofs = row_data.getIndices().size();
178 int *indices_ptr = &row_data.getIndices()[0];
179
180 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
181 iNdices.resize(nb_dofs, false);
182 noalias(iNdices) = row_data.getIndices();
184 indices_ptr = &iNdices[0];
185 }
186 frontIndices.resize(nb_dofs, false);
187 noalias(frontIndices) = row_data.getIndices();
188 VectorDofs &dofs = row_data.getFieldDofs();
189 VectorDofs::iterator dit = dofs.begin();
190 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
191 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) !=
193 iNdices[ii] = -1;
194 } else {
195 frontIndices[ii] = -1;
196 }
197 }
198 if (smootherData.frontF) {
199 CHKERR VecSetValues(smootherData.frontF, nb_dofs, &frontIndices[0],
200 &nf[0], ADD_VALUES);
201 }
202 }
203
204 CHKERR VecSetOption(getFEMethod()->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
205 PETSC_TRUE);
206 CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
207 ADD_VALUES);
208
210 }
211 };
212
215
218
220 const std::string vel_field, const std::string field_name,
223 SmootherBlockData &smoother_data,
224 const std::string
225 crack_area_tangent_constrains // = "LAMBDA_CRACK_TANGENT_CONSTRAIN"
226 )
228 data, common_data),
229 smootherData(smoother_data),
230 fieldCrackAreaTangentConstrains(crack_area_tangent_constrains) {}
231
232 ublas::vector<int> rowFrontIndices;
233
234 MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
235 EntityType col_type,
236 EntitiesFieldData::EntData &row_data,
237 EntitiesFieldData::EntData &col_data) {
239
240 int nb_row = row_data.getIndices().size();
241 int nb_col = col_data.getIndices().size();
242 int *row_indices_ptr = &row_data.getIndices()[0];
243 int *col_indices_ptr = &col_data.getIndices()[0];
244
245 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
246 rowIndices.resize(nb_row, false);
247 noalias(rowIndices) = row_data.getIndices();
249 row_indices_ptr = &rowIndices[0];
250 }
251 rowFrontIndices.resize(nb_row, false);
252 noalias(rowFrontIndices) = row_data.getIndices();
253 VectorDofs &dofs = row_data.getFieldDofs();
254 VectorDofs::iterator dit = dofs.begin();
255 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
256 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) !=
258 rowIndices[ii] = -1;
259 } else {
260 rowFrontIndices[ii] = -1;
261 }
262 }
263 }
264
265 CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr,
266 nb_col, col_indices_ptr, &k(0, 0), ADD_VALUES);
267
269
270 const auto bit_number_for_crack_area_tangent_constrain =
272 const auto bit_number_for_mesh_position =
273 getFEMethod()->getFieldBitNumber("MESH_NODE_POSITIONS");
274
275 // get tangent vector array
276 double *f_tangent_front_mesh_array;
278 &f_tangent_front_mesh_array);
279
280 auto row_dofs = getFEMethod()->getRowDofsPtr();
281
282 // iterate nodes on tet
283 for (int nn = 0; nn < 4; nn++) {
284
285 // get indices with Lagrange multiplier at node nn
286 auto dit = row_dofs->get<Unique_mi_tag>().lower_bound(
287 DofEntity::getLoFieldEntityUId(
288 bit_number_for_crack_area_tangent_constrain, getConn()[nn]));
289 auto hi_dit = row_dofs->get<Unique_mi_tag>().upper_bound(
290 DofEntity::getHiFieldEntityUId(
291 bit_number_for_crack_area_tangent_constrain, getConn()[nn]));
292
293 // continue if Lagrange are on element
294 if (std::distance(dit, hi_dit) > 0) {
295
296 // get mesh node positions at node nn
297 auto diit = row_dofs->get<Unique_mi_tag>().lower_bound(
298 DofEntity::getLoFieldEntityUId(bit_number_for_mesh_position,
299 getConn()[nn]));
300
301 auto hi_diit = row_dofs->get<Unique_mi_tag>().upper_bound(
302 DofEntity::getHiFieldEntityUId(bit_number_for_mesh_position,
303 getConn()[nn]));
304
305 // iterate over dofs on node nn
306 for (; diit != hi_diit; diit++) {
307 // iterate overt dofs in element column
308 for (int ddd = 0; ddd < nb_col; ddd++) {
309 // check consistency, node has to be at crack front
310 if (rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()] !=
311 diit->get()->getPetscGlobalDofIdx()) {
312 SETERRQ2(
313 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
314 "data inconsistency %d != %d",
315 rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()],
316 diit->get()->getPetscGlobalDofIdx());
317 }
318 // dof is not on this partition
319 if (diit->get()->getPetscLocalDofIdx() == -1) {
320 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
321 "data inconsistency");
322 }
323 double g =
324 f_tangent_front_mesh_array[diit->get()
325 ->getPetscLocalDofIdx()] *
326 k(3 * nn + diit->get()->getDofCoeffIdx(), ddd);
327 int lambda_idx = dit->get()->getPetscGlobalDofIdx();
328 CHKERR MatSetValues(getFEMethod()->snes_B, 1, &lambda_idx, 1,
329 &col_indices_ptr[ddd], &g, ADD_VALUES);
330 }
331 }
332 }
333 }
334 CHKERR VecRestoreArray(smootherData.tangentFrontF,
335 &f_tangent_front_mesh_array);
336 }
337
339 }
340 };
341};
342
343#endif //__SMOOTHER_HPP__
static PetscErrorCode ierr
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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
std::vector< MatrixDouble3by3 > sTress
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
SmootherBlockData & smootherData
Definition: Smoother.hpp:57
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: Smoother.hpp:63
MyVolumeFE(MoFEM::Interface &m_field, SmootherBlockData &smoother_data)
Definition: Smoother.hpp:59
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: Smoother.hpp:90
OpJacobianSmoother(const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, int tag, bool jacobian)
Definition: Smoother.hpp:134
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
Definition: Smoother.hpp:141
SmootherBlockData & smootherData
Definition: Smoother.hpp:216
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:219
ublas::vector< int > rowFrontIndices
Definition: Smoother.hpp:232
const std::string fieldCrackAreaTangentConstrains
Definition: Smoother.hpp:217
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:234
SmootherBlockData & smootherData
Definition: Smoother.hpp:161
ublas::vector< int > frontIndices
Definition: Smoother.hpp:171
OpRhsSmoother(const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, SmootherBlockData &smoother_data)
Definition: Smoother.hpp:163
MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: Smoother.hpp:173
MoFEMErrorCode getOptions()
Definition: Smoother.hpp:26
SmootherBlockData smootherData
Definition: Smoother.hpp:50
MyVolumeFE & getLoopFeLhs()
get lhs volume element
Definition: Smoother.hpp:124
Smoother(MoFEM::Interface &m_field)
Definition: Smoother.hpp:126
MyVolumeFE & feLhs
Definition: Smoother.hpp:123
NonlinearElasticElement::CommonData commonData
Definition: Smoother.hpp:53
boost::shared_ptr< MyVolumeFE > feRhsPtr
Definition: Smoother.hpp:118
std::map< int, NonlinearElasticElement::BlockData > setOfBlocks
Definition: Smoother.hpp:52
MyVolumeFE & feRhs
calculate right hand side for tetrahedral elements
Definition: Smoother.hpp:121
boost::shared_ptr< MyVolumeFE > feLhsPtr
Definition: Smoother.hpp:119
MyVolumeFE & getLoopFeRhs()
get rhs volume element
Definition: Smoother.hpp:122