Alamo
SutureCrack.H
Go to the documentation of this file.
1#ifndef INTEGRATOR_SUTURECRACK_H
2#define INTEGRATOR_SUTURECRACK_H
3
4#include <iostream>
5#include <fstream>
6#include <iomanip>
7
8#include "AMReX.H"
9#include "AMReX_ParallelDescriptor.H"
10#include "AMReX_ParmParse.H"
11
12#include "Integrator/Integrator.H"
13
14#include "BC/BC.H"
15#include "BC/Constant.H"
16
17#include "IC/IC.H"
18#include "IC/Laminate.H"
19#include "IC/Notch.H"
20
21#include "Operator/Operator.H"
22#include "Operator/Elastic.H"
23#include "Solver/Nonlocal/Linear.H"
24#include "Solver/Nonlocal/Newton.H"
25
26#include "Model/Solid/Solid.H"
27#include "Model/Solid/Linear/Isotropic.H"
28#include "Model/Solid/Linear/IsotropicDegradable.H"
29//#include "Model/Solid/Affine/J2PlasticDegradable.H"
30//#include "Model/Solid/Affine/CrystalPlasticDegradable.H"
31
32#include "Model/Interface/Crack/Crack.H"
33#include "Model/Interface/Crack/Constant.H"
34#include "Model/Interface/Crack/Sin.H"
35
36#include "Numeric/Stencil.H"
37#include <eigen3/Eigen/Dense>
38
39namespace Integrator
40{
41//using brittle_fracture_model_type_test = Model::Solid::Linear::Isotropic;
42using suture_fracture_model_type = Model::Solid::Linear::IsotropicDegradable;
43
44class SutureCrack : public Integrator
45{
46
47public:
48 SutureCrack()
49 {
50 IO::ParmParse pp_crack("crack");
51 pp_crack.query("modulus_scaling_max",crack.scaleModulusMax);
52 pp_crack.query("refinement_threshold",crack.refinement_threshold);
53 pp_crack.query("df_tol_rel",crack.driving_force_tolerance_rel);
54 pp_crack.query("df_tol_abs",crack.driving_force_tolerance_abs);
55
56 Model::Interface::Crack::Constant *tmpbdy = new Model::Interface::Crack::Constant();
57 pp_crack.queryclass("constant",*tmpbdy);
58 crack.cracktype = tmpbdy;
59
60 // IC for crack field
61 IO::ParmParse pp("crack.ic");
62 pp_query("type", crack.ic_type); // Crack type to use {notch}
63 if (crack.ic_type == "notch")
64 {
65 IC::Notch *tmpic = new IC::Notch(geom);
66 pp_queryclass("notch",*tmpic);
67 crack.ic = tmpic;
68 crack.is_ic = true;
69 }
70 else
71 crack.is_ic = false;
72
73 RegisterNodalFab(crack.field, 1, number_of_ghost_nodes, "crack", true);
74 RegisterNodalFab(crack.field_old, 1, number_of_ghost_nodes, "crack_old", true);
75 RegisterNodalFab(crack.driving_force, 4, number_of_ghost_nodes, "driving_force", true);
76 RegisterIntegratedVariable(&(crack.driving_force_norm),"driving_force_norm");
77
78 IO::ParmParse pp_material("material");
79 pp_material.queryclass("isotropic",material.brittlemodeltype);
80 pp_material.query("refinement_threshold",material.refinement_threshold);
81
82 IO::ParmParse pp_material_ic("material.ic");
83 pp_material_ic.query("type",material.ic_type);
84 if (material.ic_type == "laminate")
85 {
86 IC::Laminate *tmpic = new IC::Laminate(geom);
87 pp_material_ic.queryclass("laminate",*tmpic);
88
89 IO::ParmParse _temp("material.ic.laminate");
90 _temp.query("number_of_inclusions",number_of_materials);
91 number_of_materials++;
92
93 material.ic = tmpic;
94 material.is_ic = true;
95 }
96 else
97 material.is_ic = false;
98
99 IO::ParmParse pp_load("loading");
100 if (pp_load.countval("body_force")) pp_load.queryarr("body_force",loading.body_force);
101 pp_load.query("val", loading.val);
102
103 IO::ParmParse pp_elastic("elastic");
104 pp_elastic.query("df_mult", elastic.df_mult);
105 pp_elastic.queryclass("bc",elastic.brittlebc);
106
107 IO::ParmParse pp_solver("solver");
108 pp_solver.query("int", sol.interval);
109 pp_solver.query("type", sol.type);
110 pp_solver.query("max_iter", sol.max_iter);
111 pp_solver.query("max_fmg_iter", sol.max_fmg_iter);
112 pp_solver.query("verbose", sol.verbose);
113 pp_solver.query("cgverbose", sol.cgverbose);
114 pp_solver.query("tol_rel", sol.tol_rel);
115 pp_solver.query("tol_abs", sol.tol_abs);
116 pp_solver.query("cg_tol_rel", sol.cg_tol_rel);
117 pp_solver.query("cg_tol_abs", sol.cg_tol_abs);
118 pp_solver.query("use_fsmooth", sol.use_fsmooth);
119 pp_solver.query("agglomeration", sol.agglomeration);
120 pp_solver.query("consolidation", sol.consolidation);
121 pp_solver.query("bottom_solver", sol.bottom_solver);
122 pp_solver.query("linop_maxorder", sol.linop_maxorder);
123 pp_solver.query("max_coarsening_level",sol.max_coarsening_level);
124 pp_solver.query("verbose", sol.verbose);
125 pp_solver.query("cg_verbose", sol.cgverbose);
126 pp_solver.query("bottom_max_iter", sol.bottom_max_iter);
127 pp_solver.query("max_fixed_iter", sol.max_fixed_iter);
128 pp_solver.query("bottom_tol", sol.bottom_tol);
129 pp_solver.query("pre_smooth", sol.pre_smooth);
130 pp_solver.query("post_smooth", sol.post_smooth);
131
132
133 nlevels = maxLevel() + 1;
134 RegisterNodalFab(elastic.disp, AMREX_SPACEDIM, number_of_ghost_nodes, "disp", true);
135 RegisterNodalFab(elastic.rhs, AMREX_SPACEDIM, number_of_ghost_nodes, "rhs", true);
136 RegisterNodalFab(elastic.residual, AMREX_SPACEDIM, number_of_ghost_nodes, "res", true);
137 RegisterNodalFab(elastic.strain, AMREX_SPACEDIM*AMREX_SPACEDIM, number_of_ghost_nodes, "strain", true);
138 RegisterNodalFab(elastic.stress, AMREX_SPACEDIM*AMREX_SPACEDIM, number_of_ghost_nodes, "stress", true);
139 RegisterNodalFab(elastic.energy, 1, number_of_ghost_nodes, "energy", true);
140 RegisterNodalFab(elastic.energy_pristine, 1, number_of_ghost_nodes, "energy_pristine", true);
141 RegisterNodalFab(elastic.energy_pristine_old, 1, number_of_ghost_nodes, "energy_pristine_old", true);
142 // RegisterNewFab(material.modulus_field, crack.bc, 1, number_of_ghost_cells, "modulus_field", true);
143 RegisterNodalFab(material.modulus_field, number_of_materials, number_of_ghost_nodes, "modulus_field", true);
144 RegisterGeneralFab(material.brittlemodel, 1, number_of_ghost_nodes);
145
146 }
147
148protected:
149 void Initialize(int ilev) override
150 {
151 elastic.disp[ilev]->setVal(0.0);
152 elastic.strain[ilev]->setVal(0.0);
153 elastic.stress[ilev]->setVal(0.0);
154 elastic.rhs[ilev]->setVal(0.0);
155 elastic.energy[ilev]->setVal(0.0);
156 elastic.residual[ilev]->setVal(0.0);
157 elastic.energy_pristine[ilev] -> setVal(0.);
158 elastic.energy_pristine_old[ilev] -> setVal(0.);
159
160 //crack.ic->Initialize(ilev,material.modulus_field);
161 if(material.is_ic) material.ic -> Initialize(ilev,material.modulus_field);
162 else material.modulus_field[ilev]->setVal(1.0);
163
164 if(crack.is_ic)
165 {
166 crack.ic->Initialize(ilev,crack.field);
167 crack.ic->Initialize(ilev,crack.field_old);
168 }
169 else
170 {
171 crack.field[ilev]->setVal(1.0);
172 crack.field_old[ilev]->setVal(1.0);
173 }
174
175 material.brittlemodel[ilev]->setVal(material.brittlemodeltype);
176 }
177
178 void TimeStepBegin(Set::Scalar /*time*/, int iter) override
179 {
180 Util::Message(INFO,crack.driving_force_norm," ",crack.driving_force_reference," ",crack.driving_force_tolerance_rel);
181 if (crack.driving_force_norm / crack.driving_force_reference < crack.driving_force_tolerance_rel)
182 elastic.do_solve_now = true;
183 if (crack.driving_force_norm < crack.driving_force_tolerance_abs)
184 elastic.do_solve_now = true;
185
186 if (!elastic.do_solve_now) return;
187 // if(iter%sol.interval) return;
188 material.brittlemodeltype.DegradeModulus(0.0);
189
190 for (int ilev = 0; ilev < nlevels; ++ilev)
191 {
192 std::swap(elastic.energy_pristine_old[ilev], elastic.energy_pristine[ilev]);
193 material.brittlemodel[ilev]->setVal(material.brittlemodeltype);
194 Util::RealFillBoundary(*material.brittlemodel[ilev],geom[ilev]);
195 }
196
197 // Degrading the modulus
198 for (int ilev = 0; ilev < nlevels; ++ilev)
199 {
200 Util::RealFillBoundary(*crack.field[ilev],geom[ilev]);
201 Util::RealFillBoundary(*crack.field_old[ilev],geom[ilev]);
202 Util::RealFillBoundary(*material.modulus_field[ilev],geom[ilev]);
203
204 for (amrex::MFIter mfi(*elastic.disp[ilev],true); mfi.isValid(); ++mfi)
205 {
206 amrex::Box box = mfi.grownnodaltilebox();
207 amrex::Array4<const Set::Scalar> const& c_new = (*crack.field[ilev]).array(mfi);
208 amrex::Array4<const Set::Scalar> const& modbox = (*material.modulus_field[ilev]).array(mfi);
209 amrex::Array4<brittle_fracture_model_type_test> modelfab = (material.brittlemodel)[ilev]->array(mfi);
210
211 amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k){
212 Set::Scalar _temp = 0;
213 // _temp = std::min (crack.cracktype->g_phi(c_new(i,j,k,0),0), crack.cracktype->g_phi(modbox(i,j,k,0),0));
214 // );
215 // _temp = std::sqrt((crack.cracktype->g_phi(c_new(i,j,k,0),0))*(crack.cracktype->g_phi(c_new(i,j,k,0),0))
216 // + crack.cracktype->g_phi(modbox(i,j,k,0),0)*crack.cracktype->g_phi(modbox(i,j,k,0),0));
217 _temp = crack.cracktype->g_phi(c_new(i,j,k,0),0);
218 // _temp *= crack.cracktype->g_phi(modbox(i,j,k,0),0);
219 // _temp = crack.cracktype->g_phi(Numeric::Interpolate::CellToNodeAverage(c_old,i,j,k,0),0.);
220 if (std::isnan(_temp)) Util::Abort(INFO);
221 if(_temp < 0.0) Util::Abort(INFO);//_temp = 0.;
222 if(_temp > 1.0) Util::Abort(INFO);//_temp = 1.0;
223 _temp = crack.scaleModulusMax + _temp * (1. - crack.scaleModulusMax);
224 //modelfab(i,j,k,0).DegradeModulus(std::min(1.-_temp,1.-crack.scaleModulusMax));
225 modelfab(i,j,k,0).DegradeModulus(1-_temp);
226 });
227 }
228 Util::RealFillBoundary(*material.brittlemodel[ilev],geom[ilev]);
229 }
230
231 // Body force
232 for (int ilev = 0; ilev < nlevels; ++ilev)
233 {
234 const Real* DX = geom[ilev].CellSize();
235 Set::Scalar volume = AMREX_D_TERM(DX[0],*DX[1],*DX[2]);
236
237 AMREX_D_TERM(elastic.rhs[ilev]->setVal(loading.body_force(0)*volume,0,1);,
238 elastic.rhs[ilev]->setVal(loading.body_force(1)*volume,1,1);,
239 elastic.rhs[ilev]->setVal(loading.body_force(2)*volume,2,1););
240 }
241
242 elastic.brittlebc.Init(elastic.rhs,geom);
243
244 Operator::Elastic<brittle_fracture_model_type_test::sym> op_b;
245 LPInfo info;
246 info.setAgglomeration(sol.agglomeration);
247 info.setConsolidation(sol.consolidation);
248 info.setMaxCoarseningLevel(sol.max_coarsening_level);
249
250 for (int ilev = 0; ilev < nlevels; ilev++) if (elastic.disp[ilev]->contains_nan()) Util::Warning(INFO);
251
252 op_b.define(geom, grids, dmap, info);
253 op_b.setMaxOrder(sol.linop_maxorder);
254 op_b.SetBC(&elastic.brittlebc);
255
256 {
257 Solver::Nonlocal::Newton<brittle_fracture_model_type_test> solver(op_b);
258 solver.setMaxIter(sol.max_iter);
259 solver.setMaxFmgIter(sol.max_fmg_iter);
260 solver.setFixedIter(sol.max_fixed_iter);
261 solver.setVerbose(sol.verbose);
262 //solver.setCGVerbose(sol.cgverbose);
263 solver.setBottomMaxIter(sol.bottom_max_iter);
264 solver.setBottomTolerance(sol.cg_tol_rel) ;
265 solver.setBottomToleranceAbs(sol.cg_tol_abs) ;
266 solver.setPreSmooth(sol.pre_smooth);
267 solver.setPostSmooth(sol.post_smooth);
268 if (sol.bottom_solver == "cg") solver.setBottomSolver(MLMG::BottomSolver::cg);
269 else if (sol.bottom_solver == "bicgstab") solver.setBottomSolver(MLMG::BottomSolver::bicgstab);
270 solver.solve(elastic.disp, elastic.rhs, material.brittlemodel, sol.tol_rel, sol.tol_abs);
271 solver.compResidual(elastic.residual,elastic.disp,elastic.rhs,material.brittlemodel);
272
273 //if (iter == 0)
274 //{
275 // solver.solve(elastic.disp, elastic.rhs, material.brittlemodel, sol.tol_rel, sol.tol_abs);
276 // solver.compResidual(elastic.residual,elastic.disp,elastic.rhs,material.brittlemodel);
277 //}
278 }
279
280 for (int ilev = 0; ilev < nlevels; ilev++)
281 {
282 op_b.Strain(ilev,*elastic.strain[ilev],*elastic.disp[ilev]);
283 op_b.Stress(ilev,*elastic.stress[ilev],*elastic.disp[ilev]);
284 op_b.Energy(ilev,*elastic.energy[ilev],*elastic.disp[ilev]);
285
286 Util::RealFillBoundary(*elastic.strain[ilev],geom[ilev]);
287 Util::RealFillBoundary(*elastic.energy_pristine_old[ilev],geom[ilev]);
288 Util::RealFillBoundary(*material.modulus_field[ilev],geom[ilev]);
289 elastic.energy_pristine[ilev]->setVal(0.0);
290
291 for (amrex::MFIter mfi(*elastic.strain[ilev],true); mfi.isValid(); ++mfi)
292 {
293 const amrex::Box& box = mfi.grownnodaltilebox();
294 amrex::Array4<const Set::Scalar> const& strain_box = (*elastic.strain[ilev]).array(mfi);
295 // amrex::Array4<const Set::Scalar> const& modbox = (*material.modulus_field[ilev]).array(mfi);
296 amrex::Array4<Set::Scalar> const& energy_box = (*elastic.energy_pristine[ilev]).array(mfi);
297 amrex::Array4<Set::Scalar> const& energy_box_old = (*elastic.energy_pristine_old[ilev]).array(mfi);
298
299 amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k){
300 Set::Matrix eps = Numeric::FieldToMatrix(strain_box,i,j,k);
301 Eigen::SelfAdjointEigenSolver<Set::Matrix> eigensolver(eps);
302 Set::Vector eValues = eigensolver.eigenvalues();
303 Set::Matrix eVectors = eigensolver.eigenvectors();
304
305 Set::Matrix epsp = Set::Matrix::Zero();
306 Set::Matrix epsn = Set::Matrix::Zero();
307
308 for (int n = 0; n < AMREX_SPACEDIM; n++)
309 {
310 if(eValues(n) > 0.0) epsp += eValues(n)*(eVectors.col(n)*eVectors.col(n).transpose());
311 else epsn += eValues(n)*(eVectors.col(n)*eVectors.col(n).transpose());
312 }
313 // Set::Scalar _temp = crack.cracktype->g_phi(modbox(i,j,k,0),0.);
314 // if(std::isnan(_temp)) Util::Abort(INFO, "Nans in temp. modbox(", i,",",j,",",k,") = ", modbox(i,j,k,0));
315 // if (_temp < 0.0) _temp = 0.0;
316 // if (_temp > 1.0) _temp = 1.0;
317 // _temp = crack.scaleModulusMax + _temp * (1. - crack.scaleModulusMax);
318 // material.brittlemodeltype.DegradeModulus(1.- _temp);
319 energy_box(i,j,k,0) = material.brittlemodeltype.W(epsp);
320 energy_box(i,j,k,0) = energy_box(i,j,k,0) > energy_box_old(i,j,k,0) ? energy_box(i,j,k,0) : energy_box_old(i,j,k,0);
321 });
322 }
323 Util::RealFillBoundary(*elastic.energy_pristine[ilev],geom[ilev]);
324 }
325
326 integrate_variables_before_advance = false;
327 integrate_variables_after_advance = true;
328 }
329
330 void Advance(int lev, Set::Scalar /*time*/, Set::Scalar dt) override
331 {
332 std::swap(crack.field_old[lev], crack.field[lev]);
333
334 // crack.field_old[lev]->FillBoundary();
335 // material.modulus_field[lev]->FillBoundary();
336 Util::RealFillBoundary(*crack.field_old[lev],geom[lev]);
337 Util::RealFillBoundary(*crack.field[lev],geom[lev]);
338 Util::RealFillBoundary(*material.modulus_field[lev],geom[lev]);
339 Util::RealFillBoundary(*elastic.energy_pristine[lev],geom[lev]);
340
341 const Set::Scalar* DX = geom[lev].CellSize();
342 amrex::Box domain(geom[lev].Domain());
343 // amrex::Box domain(geom.Domain();
344 domain.convert(amrex::IntVect::TheNodeVector());
345 const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
346
347 for ( amrex::MFIter mfi(*crack.field[lev],true); mfi.isValid(); ++mfi )
348 {
349 const amrex::Box& bx = mfi.validbox();
350 amrex::Array4<const Set::Scalar> const& c_old = (*crack.field_old[lev]).array(mfi);
351 amrex::Array4<Set::Scalar> const& df = (*crack.driving_force[lev]).array(mfi);
352 amrex::Array4<Set::Scalar> const& c_new = (*crack.field[lev]).array(mfi);
353 amrex::Array4<const Set::Scalar> const& energy_box = (*elastic.energy_pristine[lev]).array(mfi);
354 amrex::Array4<const Set::Scalar> const& modbox = (*material.modulus_field[lev]).array(mfi);
355
356 amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k){
357
358#if AMREX_SPACEDIM !=2
359 Util::Abort(INFO, "This doesn't work with 1D or 3D yet");
360#endif
361 if (i == lo.x && j == lo.y) c_new(i,j,k,0) = c_new(i+1,j+1,k,0);
362 else if (i == lo.x && j == hi.y) c_new(i,j,k,0) = c_new(i+1,j-1,k,0);
363 else if (i == hi.x && j == lo.y) c_new(i,j,k,0) = c_new(i-1,j+1,k,0);
364 else if (i == hi.x && j == hi.y) c_new(i,j,k,0) = c_new(i-1,j-1,k,0);
365 else if (i == lo.x) c_new(i,j,k,0) = c_new(i+1,j,k,0);
366 else if (j == lo.y) c_new(i,j,k,0) = c_new(i,j+1,k,0);
367 else if (i == hi.x) c_new(i,j,k,0) = c_new(i-1,j,k,0);
368 else if (j == hi.y) c_new(i,j,k,0) = c_new(i,j-1,k,0);
369 else
370 {
371 Set::Scalar rhs = 0.0;
372 //Set::Vector Dc = Numeric::Gradient(c_old, i, j, k, 0, DX);
373
374 Set::Scalar en_cell = energy_box(i,j,k,0);
375
376 if (std::isnan(en_cell)) Util::Abort(INFO, "Nans detected in en_cell. energy_box(i,j,k,0) = ", energy_box(i,j,k,0));
377 if (std::isinf(c_old(i,j,k,0))) Util::Abort(INFO, "Infs detected in c_old");
378 //if (c_old(i,j,k,0) > 1) Util::Abort(INFO, "Very large values of c_old at (",i,",",j,",",k,") = ", c_old(i,j,k,0));
379
380 df(i,j,k,0) = crack.cracktype->Dg_phi(c_old(i,j,k,0),0.0)*en_cell*elastic.df_mult;
381 rhs += crack.cracktype->Dg_phi(c_old(i,j,k,0),0.0)*en_cell*elastic.df_mult;
382
383 //if(rhs > 1.e10) Util::Abort(INFO, "Very large values of rhs at (",i,",",j,",",k,") = ", rhs);
384
385 Set::Matrix DDc = Numeric::Hessian(c_old, i, j, k, 0, DX);
386 Set::Scalar laplacian = DDc.trace();
387
388 Set::Scalar _temp = 1.0;
389 for (int m = 0; m < number_of_materials; m++)
390 _temp *= modbox(i,j,k,m);
391
392 Set::Scalar Gc = (1.0 - _temp)*crack.cracktype->Gc(0.0);
393
394 df(i,j,k,1) = Gc*crack.cracktype->Dw_phi(c_old(i,j,k,0),0.0)/(4.0*crack.cracktype->Zeta(0.0))*crack.mult_df_Gc;
395 df(i,j,k,2) = 2.0*crack.cracktype->Zeta(0.0)*Gc*laplacian*crack.mult_df_lap;
396
397 rhs += Gc*crack.cracktype->Dw_phi(c_old(i,j,k,0),0.)/(4.0*crack.cracktype->Zeta(0.0))*crack.mult_df_Gc;
398 // if(rhs > 1.e10) Util::Abort(INFO, "Very large values of rhs at (",i,",",j,",",k,") = ", rhs, ". c_old(i,j,k,0) = ",c_old(i,j,k,0));
399
400 rhs -= 2.0*crack.cracktype->Zeta(0.0)*Gc*laplacian*crack.mult_df_lap;
401 // if(rhs > 1.e10) Util::Abort(INFO, "Very large values of rhs at (",i,",",j,",",k,") = ", rhs, ". c_old(i,j,k,0) = ",c_old(i,j,k,0), ". laplacian = ", laplacian);
402
403 // rhs *= _temp;
404 // rhs *= crack.cracktype->g_phi(_temp,0.0);
405 // if(rhs > 1.e10) Util::Abort(INFO, "Very large values of rhs at (",i,",",j,",",k,") = ", rhs, ". c_old(i,j,k,0) = ",c_old(i,j,k,0));
406
407 df(i,j,k,3) = std::max(0.,rhs - crack.cracktype->DrivingForceThreshold(c_old(i,j,k,0)));
408 if(std::isnan(rhs)) Util::Abort(INFO, "Dwphi = ", crack.cracktype->Dw_phi(c_old(i,j,k,0),0.0),". c_old(i,j,k,0) = ",c_old(i,j,k,0));
409 c_new(i,j,k,0) = c_old(i,j,k,0) - dt*std::max(0., rhs - crack.cracktype->DrivingForceThreshold(c_old(i,j,k,0)))*crack.cracktype->Mobility(c_old(i,j,k,0));
410
411 if (c_new (i,j,k,0) < 0.0) c_new(i,j,k,0) = 0.0;
412 if (c_new (i,j,k,0) > 1.0) c_new(i,j,k,0) = 1.0;
413 }
414 });
415 }
416 // crack.field[lev]->FillBoundary();
417 Util::RealFillBoundary(*crack.field[lev],geom[lev]);
418 }
419
420 void TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, amrex::Real /*time*/, int /*ngrow*/) override
421 {
422 const amrex::Real *DX = geom[lev].CellSize();
423 const Set::Vector dx(DX);
424 const Set::Scalar dxnorm = dx.lpNorm<2>();
425 amrex::Box domain(geom[lev].Domain());
426 domain.convert(amrex::IntVect::TheNodeVector());
427
428 for (amrex::MFIter mfi(*crack.field[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
429 {
430 const amrex::Box bx = mfi.tilebox() & domain.grow(-1);
431 amrex::Array4<char> const &tags = a_tags.array(mfi);
432 amrex::Array4<const Set::Scalar> const &c_new = (*crack.field[lev]).array(mfi);
433 amrex::Array4<const Set::Scalar> const &modbox = (*material.modulus_field[lev]).array(mfi);
434
435 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
436 Set::Vector grad = Numeric::Gradient(c_new, i, j, k, 0, DX);
437 Set::Vector grad2 = Numeric::Gradient(modbox, i, j, k, 0, DX);
438 if (dxnorm * grad.lpNorm<2>() >= crack.refinement_threshold || dxnorm * grad2.lpNorm<2>() >= material.refinement_threshold)
439 tags(i, j, k) = amrex::TagBox::SET;
440 });
441 }
442 }
443
444 void Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,const amrex::MFIter &mfi, const amrex::Box &box)
445 {
446 const amrex::Real* DX = geom[amrlev].CellSize();
447 const Set::Scalar DV = AMREX_D_TERM(DX[0],*DX[1],*DX[2]);
448
449 amrex::Array4<const Set::Scalar> const &df = (*crack.driving_force[amrlev]).array(mfi);
450
451 amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
452 {
453 crack.driving_force_norm += df(i,j,k,3) * DV;
454 });
455 }
456
457 void TimeStepComplete(amrex::Real /*time*/,int /*iter*/)
458 {
459 if (elastic.do_solve_now)
460 crack.driving_force_reference = crack.driving_force_norm;
461
462 elastic.do_solve_now = false;
463 }
464
465private:
466
467 int number_of_ghost_nodes = 2; ///< Number of ghost nodes
468 int number_of_materials = 1;
469 int nlevels;
470
471 struct{
472 Set::Field<Set::Scalar> disp; ///< displacement field
473 Set::Field<Set::Scalar> strain; ///< total strain field (gradient of displacement)
474 Set::Field<Set::Scalar> stress; ///< stress field
475 Set::Field<Set::Scalar> rhs; ///< rhs fab for elastic solution
476 Set::Field<Set::Scalar> residual; ///< residual field for solver
477 Set::Field<Set::Scalar> energy; ///< total elastic energy
478 Set::Field<Set::Scalar> energy_pristine; ///< energy of the prisitne material as if no crack is present
479 Set::Field<Set::Scalar> energy_pristine_old; ///< energy of the pristine material for previous time step.
480
481 BC::Operator::Elastic::Constant brittlebc; ///< elastic BC if using brittle fracture
482 Set::Scalar df_mult = 1.0; ///< mulitplier for elastic driving force.
483 bool do_solve_now = false;
484 } elastic;
485
486 struct{
487 Set::Field<Set::Scalar> field; ///< crack field at current time step
488 Set::Field<Set::Scalar> field_old; ///< crack field at previous time step
489 Set::Field<Set::Scalar> driving_force; ///< crack driving forces.
490 Set::Scalar driving_force_reference = 1.0;
491 Set::Scalar driving_force_norm = 0.0;
492 Set::Scalar driving_force_tolerance_rel = 1E-4;
493 Set::Scalar driving_force_tolerance_abs = 0.0;
494
495 Model::Interface::Crack::Crack *cracktype; ///< type of crack. See Crack/Constant or Crack/Sin
496 std::string ic_type; ///< crack IC type. See IC/Notch and IC/Ellipsoid
497 IC::IC *ic; ///< crack IC. See IC/Notch and IC/Ellipsoid
498 bool is_ic = false;
499
500 Set::Scalar scaleModulusMax = 0.02; ///< material modulus ratio inside crack (default = 0.02).
501 Set::Scalar refinement_threshold = 0.001;
502
503 Set::Scalar mult_df_Gc = 1.0; ///< Multiplier for Gc/zeta term
504 Set::Scalar mult_df_lap = 1.0; ///< Multiplier for the laplacian term
505 } crack;
506
507 struct{
508 suture_fracture_model_type brittlemodeltype;
509 Set::Field<Set::Scalar> modulus_field;
510 Set::Field<suture_fracture_model_type> brittlemodel;
511 std::string input_material = "isotropic";
512 Set::Scalar refinement_threshold = 0.1;
513
514 IC::IC *ic;
515 std::string ic_type;
516 bool is_ic = false;
517 } material;
518
519 struct{
520 Set::Scalar bottom_tol = 1.0E-6;
521 int interval = 100;
522 std::string type = "single";
523 int max_iter = 1000;
524 int max_fmg_iter = 0;
525 int bottom_max_iter = 200;
526 int max_fixed_iter = 500;
527 int verbose = 3;
528 int cgverbose = 3;
529 Set::Scalar tol_rel = 1.0E-6;
530 Set::Scalar tol_abs = 1.0E-6;
531 Set::Scalar cg_tol_rel = 1.0E-6;
532 Set::Scalar cg_tol_abs = 1.0E-6;
533 Set::Scalar tstart = 0.0;
534 Set::Scalar tend = 0.0;
535 std::string bottom_solver = "bicgstab";
536 int linop_maxorder = 2;
537 bool use_fsmooth = false;
538 int max_coarsening_level = 0;
539 bool agglomeration = true;
540 bool consolidation = false;
541 int pre_smooth = 2;
542 int post_smooth = 2;
543 } sol;
544
545 struct{
546 Set::Vector body_force = Set::Vector::Zero();
547 Set::Scalar val = 0.;
548 } loading;
549
550};
551}
552#endif
pp_query
#define pp_query(...)
Definition ParmParse.H:106
pp_queryclass
#define pp_queryclass(...)
Definition ParmParse.H:107
INFO
#define INFO
Definition Util.H:20
IC::IC
Pure abstract IC object from which all other IC objects inherit.
Definition IC.H:23
IC::Laminate
Initialize Laminates in a matrix.
Definition Laminate.H:16
IO::ParmParse::queryclass
void queryclass(std::string name, T *value, std::string file="", std::string func="", int line=-1)
Definition ParmParse.H:604
IO::ParmParse::queryarr
int queryarr(std::string name, std::vector< T > &value, std::string, std::string, int)
Definition ParmParse.H:336
Integrator::Integrator::box
std::vector< amrex::Box > box
Definition Integrator.H:446
Integrator::Integrator::RegisterNodalFab
void RegisterNodalFab(Set::Field< Set::Scalar > &new_fab, int ncomp, int nghost, std::string name, bool writeout, std::vector< std::string > suffix={})
Add a new node-based scalar field.
Definition Integrator.cpp:345
Integrator::Integrator::dt
amrex::Vector< amrex::Real > dt
Timesteps for each level of refinement.
Definition Integrator.H:381
Integrator::Integrator::integrate_variables_before_advance
bool integrate_variables_before_advance
Definition Integrator.H:392
Integrator::Integrator::RegisterGeneralFab
void RegisterGeneralFab(Set::Field< T > &new_fab, int ncomp, int nghost, bool evolving=true)
Add a templated nodal field.
Integrator::Integrator::RegisterIntegratedVariable
void RegisterIntegratedVariable(Set::Scalar *integrated_variable, std::string name, bool extensive=true)
Register a variable to be integrated over the spatial domain using the Integrate function.
Definition Integrator.cpp:355
Integrator::Integrator::integrate_variables_after_advance
bool integrate_variables_after_advance
Definition Integrator.H:393
Integrator::SutureCrack::TimeStepComplete
void TimeStepComplete(amrex::Real, int)
Definition SutureCrack.H:457
Integrator::SutureCrack::brittlemodeltype
suture_fracture_model_type brittlemodeltype
Definition SutureCrack.H:508
Integrator::SutureCrack::field
Set::Field< Set::Scalar > field
crack field at current time step
Definition SutureCrack.H:487
Integrator::SutureCrack::ic
IC::IC * ic
crack IC. See IC/Notch and IC/Ellipsoid
Definition SutureCrack.H:497
Integrator::SutureCrack::elastic
struct Integrator::SutureCrack::@26 elastic
Integrator::SutureCrack::field_old
Set::Field< Set::Scalar > field_old
crack field at previous time step
Definition SutureCrack.H:488
Integrator::SutureCrack::driving_force_norm
Set::Scalar driving_force_norm
Definition SutureCrack.H:491
Integrator::SutureCrack::brittlemodel
Set::Field< suture_fracture_model_type > brittlemodel
Definition SutureCrack.H:510
Integrator::SutureCrack::rhs
Set::Field< Set::Scalar > rhs
rhs fab for elastic solution
Definition SutureCrack.H:475
Integrator::SutureCrack::driving_force_tolerance_rel
Set::Scalar driving_force_tolerance_rel
Definition SutureCrack.H:492
Integrator::SutureCrack::driving_force_reference
Set::Scalar driving_force_reference
Definition SutureCrack.H:490
Integrator::SutureCrack::modulus_field
Set::Field< Set::Scalar > modulus_field
Definition SutureCrack.H:509
Integrator::SutureCrack::crack
struct Integrator::SutureCrack::@27 crack
Integrator::SutureCrack::scaleModulusMax
Set::Scalar scaleModulusMax
material modulus ratio inside crack (default = 0.02).
Definition SutureCrack.H:500
Integrator::SutureCrack::mult_df_Gc
Set::Scalar mult_df_Gc
Multiplier for Gc/zeta term.
Definition SutureCrack.H:503
Integrator::SutureCrack::material
struct Integrator::SutureCrack::@28 material
Integrator::SutureCrack::TimeStepBegin
void TimeStepBegin(Set::Scalar, int iter) override
Definition SutureCrack.H:178
Integrator::SutureCrack::TagCellsForRefinement
void TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, amrex::Real, int) override
Definition SutureCrack.H:420
Integrator::SutureCrack::energy
Set::Field< Set::Scalar > energy
total elastic energy
Definition SutureCrack.H:477
Integrator::SutureCrack::mult_df_lap
Set::Scalar mult_df_lap
Multiplier for the laplacian term.
Definition SutureCrack.H:504
Integrator::SutureCrack::sol
struct Integrator::SutureCrack::@29 sol
Integrator::SutureCrack::residual
Set::Field< Set::Scalar > residual
residual field for solver
Definition SutureCrack.H:476
Integrator::SutureCrack::brittlebc
BC::Operator::Elastic::Constant brittlebc
elastic BC if using brittle fracture
Definition SutureCrack.H:481
Integrator::SutureCrack::energy_pristine
Set::Field< Set::Scalar > energy_pristine
energy of the prisitne material as if no crack is present
Definition SutureCrack.H:478
Integrator::SutureCrack::Integrate
void Integrate(int amrlev, Set::Scalar, int, const amrex::MFIter &mfi, const amrex::Box &box)
Definition SutureCrack.H:444
Integrator::SutureCrack::refinement_threshold
Set::Scalar refinement_threshold
Definition SutureCrack.H:501
Integrator::SutureCrack::cracktype
Model::Interface::Crack::Crack * cracktype
type of crack. See Crack/Constant or Crack/Sin
Definition SutureCrack.H:495
Integrator::SutureCrack::driving_force
Set::Field< Set::Scalar > driving_force
crack driving forces.
Definition SutureCrack.H:489
Integrator::SutureCrack::Initialize
void Initialize(int ilev) override
Definition SutureCrack.H:149
Integrator::SutureCrack::driving_force_tolerance_abs
Set::Scalar driving_force_tolerance_abs
Definition SutureCrack.H:493
Integrator::SutureCrack::number_of_ghost_nodes
int number_of_ghost_nodes
Number of ghost nodes.
Definition SutureCrack.H:467
Integrator::SutureCrack::stress
Set::Field< Set::Scalar > stress
stress field
Definition SutureCrack.H:474
Integrator::SutureCrack::loading
struct Integrator::SutureCrack::@30 loading
Integrator::SutureCrack::strain
Set::Field< Set::Scalar > strain
total strain field (gradient of displacement)
Definition SutureCrack.H:473
Integrator::SutureCrack::disp
Set::Field< Set::Scalar > disp
displacement field
Definition SutureCrack.H:472
Integrator::SutureCrack::Advance
void Advance(int lev, Set::Scalar, Set::Scalar dt) override
Definition SutureCrack.H:330
Integrator::SutureCrack::energy_pristine_old
Set::Field< Set::Scalar > energy_pristine_old
energy of the pristine material for previous time step.
Definition SutureCrack.H:479
Integrator::SutureCrack::ic_type
std::string ic_type
crack IC type. See IC/Notch and IC/Ellipsoid
Definition SutureCrack.H:496
Integrator::SutureCrack::df_mult
Set::Scalar df_mult
mulitplier for elastic driving force.
Definition SutureCrack.H:482
Operator::Elastic::define
void define(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const LPInfo &a_info=LPInfo(), const Vector< FabFactory< FArrayBox > const * > &a_factory={})
Definition Elastic.cpp:26
Operator::Elastic::SetBC
void SetBC(::BC::Operator::Elastic::Elastic *a_bc)
The different types of Boundary Condtiions are listed in the BC::Operator::Elastic documentation.
Definition Elastic.H:61
Operator::Elastic::Stress
void Stress(int amrlev, amrex::MultiFab &sigmafab, const amrex::MultiFab &ufab, bool voigt=false, bool a_homogeneous=false)
Compute stress given the displacement field by.
Definition Elastic.cpp:464
Operator::Elastic::Strain
void Strain(int amrlev, amrex::MultiFab &epsfab, const amrex::MultiFab &ufab, bool voigt=false) const
Compute strain given the displacement field by.
Definition Elastic.cpp:403
Operator::Elastic::Energy
void Energy(int amrlev, amrex::MultiFab &energy, const amrex::MultiFab &u, bool a_homogeneous=false)
Compute energy density given the displacement field by.
Definition Elastic.cpp:529
Solver::Nonlocal::Linear::setPostSmooth
void setPostSmooth(const int a_post_smooth)
Definition Linear.H:185
Solver::Nonlocal::Linear::setBottomMaxIter
void setBottomMaxIter(const int a_bottom_max_iter)
Definition Linear.H:180
Solver::Nonlocal::Linear::setVerbose
void setVerbose(const int a_verbose)
Definition Linear.H:183
Solver::Nonlocal::Linear::setMaxFmgIter
void setMaxFmgIter(const int a_max_fmg_iter)
Definition Linear.H:181
Solver::Nonlocal::Linear::setFixedIter
void setFixedIter(const int a_fixed_iter)
Definition Linear.H:182
Solver::Nonlocal::Linear::setPreSmooth
void setPreSmooth(const int a_pre_smooth)
Definition Linear.H:184
Solver::Nonlocal::Linear::setMaxIter
void setMaxIter(const int a_max_iter)
Definition Linear.H:179
Solver::Nonlocal::Newton::compResidual
void compResidual(Set::Field< Set::Scalar > &a_res_mf, Set::Field< Set::Scalar > &a_u_mf, Set::Field< Set::Scalar > &a_b_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:442
Solver::Nonlocal::Newton::solve
Set::Scalar solve(const Set::Field< Set::Vector > &a_u_mf, const Set::Field< Set::Vector > &a_b_mf, Set::Field< T > &a_model_mf, Real a_tol_rel, Real a_tol_abs, const char *checkpoint_file=nullptr)
Definition Newton.H:275
Integrator
Collection of numerical integrator objects.
Definition AllenCahn.H:41
Integrator::suture_fracture_model_type
Model::Solid::Linear::IsotropicDegradable suture_fracture_model_type
Definition SutureCrack.H:42
Numeric::FieldToMatrix
AMREX_FORCE_INLINE Set::Matrix FieldToMatrix(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k)
Definition Stencil.H:1064
Numeric::Hessian
AMREX_FORCE_INLINE Set::Matrix Hessian(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition Stencil.H:946
Numeric::Gradient
AMREX_FORCE_INLINE Set::Vector Gradient(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition Stencil.H:672
Set::Scalar
amrex::Real Scalar
Definition Base.H:19
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:20
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition Base.H:23
Util::RealFillBoundary
AMREX_FORCE_INLINE void RealFillBoundary(amrex::FabArray< amrex::BaseFab< T > > &a_mf, const amrex::Geometry &)
Definition Util.H:322
Util::Abort
void Abort(const char *msg)
Definition Util.cpp:170
Util::Warning
void Warning(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:181
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:141
Generated by doxygen 1.9.8