Alamo
Fracture.H
Go to the documentation of this file.
1#ifndef INTEGRATOR_FRACTURE_H
2#define INTEGRATOR_FRACTURE_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/PerturbedInterface.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{
41using brittle_fracture_model_type_test = Model::Solid::Linear::Isotropic;
42//using brittle_fracture_model_type_test = Model::Solid::Linear::IsotropicDegradable;
43
44class Fracture : public Integrator
45{
46
47public:
48 static constexpr const char* name = "fracture";
49
50 Fracture()
51 {
52 IO::ParmParse pp_crack("crack");
53 pp_crack.query("modulus_scaling_max",crack.scaleModulusMax);
54 pp_crack.query("refinement_threshold",crack.refinement_threshold);
55
56
57 pp_crack.queryclass("constant",crack.cracktype);
58
59 IO::ParmParse pp_crack_df("crack.df");
60 pp_crack_df.query("mult_Gc", crack.mult_df_Gc);
61 pp_crack_df.query("mult_Lap", crack.mult_df_lap);
62 pp_crack_df.query("beta", crack.beta);
63 pp_crack_df.query("tol_rel",crack.driving_force_tolerance_rel);
64 pp_crack_df.query("tol_abs",crack.driving_force_tolerance_abs);
65
66 // IC for crack field
67 IO::ParmParse pp("crack.ic");
68 pp_queryarr("type", crack.ic_type); // Type of crack {notch,ellipsoid}
69
70 if (crack.ic_type.size() == 0) crack.is_ic = false;
71 else
72 {
73 crack.is_ic = true;
74 for (int i = 0; i < crack.ic_type.size(); i++)
75 {
76 if(crack.ic_type[i] == "notch")
77 {
78 Util::Abort(INFO,"Not Supported");
79 // IC::Notch *tmpic = new IC::Notch(geom);
80 // pp_queryclass("notch",*tmpic);
81 // crack.ic.push_back(tmpic);
82 }
83 else if (crack.ic_type[i] == "ellipsoid")
84 {
85 Util::Abort(INFO,"Not Supported");
86 // IC::Ellipsoid *tmpic = new IC::Ellipsoid(geom);
87 // pp_queryclass("ellipsoid",*tmpic);
88 // crack.ic.push_back(tmpic);
89 }
90 else
91 Util::Abort(INFO, "This IC hasn't been implemented yet");
92 }
93 }
94
95 RegisterNodalFab(crack.c_mf, 1, number_of_ghost_nodes, "crack", true);
96 RegisterNodalFab(crack.c_old_mf, 1, number_of_ghost_nodes, "crack_old", false);
97 RegisterNodalFab(crack.driving_force_mf, 5, number_of_ghost_nodes, "driving_force", false);
98
99 RegisterIntegratedVariable(&(crack.driving_force_norm),"driving_force_norm");
100 RegisterIntegratedVariable(&(crack.int_crack),"int_c");
101 RegisterIntegratedVariable(&(elastic.int_energy),"int_W");
102
103 IO::ParmParse pp_material("material");
104 pp_material.query("refinement_threshold",material.refinement_threshold);
105
106 IO::ParmParse pp_material_ic("material.ic");
107 pp_material_ic.query("type",material.ic_type);
108 if (material.ic_type == "laminate")
109 {
110 IC::Laminate *tmpic = new IC::Laminate(geom);
111 pp_material_ic.queryclass("laminate",*tmpic);
112
113 IO::ParmParse _temp("material.ic.laminate");
114 _temp.query("number_of_inclusions",number_of_materials);
115 number_of_materials++;
116
117 material.ic = tmpic;
118 material.is_ic = true;
119 }
120 else if (material.ic_type == "ellipse")
121 {
122 IC::Ellipse *tmpic = new IC::Ellipse(geom);
123 pp_material_ic.queryclass("ellipse",*tmpic);
124 number_of_materials = 2;
125
126 material.ic = tmpic;
127 material.is_ic = true;
128 }
129 else if (material.ic_type == "perturbed_interface")
130 {
131 IC::PerturbedInterface *tmpic = new IC::PerturbedInterface(geom);
132 pp_material_ic.queryclass("perturbed_interface", *tmpic);
133 number_of_materials = 2;
134
135 material.ic = tmpic;
136 material.is_ic = true;
137 }
138 else
139 material.is_ic = false;
140
141 for (int p = 0; p < number_of_materials; p++)
142 {
143 std::string pp_string = "material.model"+std::to_string(p+1);
144 Util::Message(INFO, pp_string);
145 IO::ParmParse pp_model(pp_string);
146
147 if(pp_model.contains("type"))
148 {
149 std::string model_type_str;
150 pp_model.query("type", model_type_str);
151 if(model_type_str == "isotropic")
152 {
153 brittle_fracture_model_type_test *tmpmodel = new brittle_fracture_model_type_test();
154 pp_model.queryclass("isotropic",*tmpmodel);
155 material.brittlemodeltype.push_back(*tmpmodel);
156 }
157 }
158 else
159 {
160 if (p == 0) Util::Abort(INFO, "Need material information for at least one material");
161 material.brittlemodeltype.push_back(material.brittlemodeltype[p-1]);
162 }
163 }
164
165 IO::ParmParse pp_load("loading");
166 if (pp_load.countval("body_force")) pp_load.queryarr("body_force",loading.body_force);
167 pp_load.query("val", loading.val);
168
169 IO::ParmParse pp_elastic("elastic");
170 pp_elastic.query("df_mult", elastic.df_mult);
171 pp_elastic.queryclass("bc",elastic.brittlebc);
172 pp_elastic.query("int",elastic.interval);
173 pp_elastic.query("omega",elastic.omega);
174
175
176 nlevels = maxLevel() + 1;
177 RegisterNodalFab(elastic.disp_mf, AMREX_SPACEDIM, number_of_ghost_nodes, "disp", true);
178 RegisterNodalFab(elastic.rhs_mf, AMREX_SPACEDIM, number_of_ghost_nodes, "rhs", false);
179 RegisterNodalFab(elastic.residual_mf, AMREX_SPACEDIM, number_of_ghost_nodes, "res", false);
180 RegisterNodalFab(elastic.strain_mf, AMREX_SPACEDIM*AMREX_SPACEDIM, number_of_ghost_nodes, "strain", false);
181 RegisterNodalFab(elastic.stress_mf, AMREX_SPACEDIM*AMREX_SPACEDIM, number_of_ghost_nodes, "stress", true);
182 RegisterNodalFab(elastic.energy_mf, 1, number_of_ghost_nodes, "energy", false);
183 RegisterNodalFab(elastic.energy_pristine_mf, 1, number_of_ghost_nodes, "energy_pristine", false);
184 RegisterNodalFab(elastic.energy_pristine_old_mf, 1, number_of_ghost_nodes, "energy_pristine_old", false);
185 RegisterNodalFab(material.material_mf, number_of_materials, number_of_ghost_nodes, "materials", true);
186 RegisterGeneralFab(material.model_mf, 1, number_of_ghost_nodes);
187
188 }
189
190protected:
191 void Initialize(int ilev) override
192 {
193 elastic.disp_mf[ilev]->setVal(0.0);
194 elastic.strain_mf[ilev]->setVal(0.0);
195 elastic.stress_mf[ilev]->setVal(0.0);
196 elastic.rhs_mf[ilev]->setVal(0.0);
197 elastic.energy_mf[ilev]->setVal(0.0);
198 elastic.residual_mf[ilev]->setVal(0.0);
199 elastic.energy_pristine_mf[ilev] -> setVal(0.);
200 elastic.energy_pristine_old_mf[ilev] -> setVal(0.);
201
202 if(crack.is_ic)
203 {
204 crack.c_mf[ilev]->setVal(1.0);
205 crack.c_old_mf[ilev]->setVal(1.0);
206
207 for (int i = 0; i < crack.ic.size(); i++)
208 {
209 crack.ic[i]->Add(ilev,crack.c_mf);
210 crack.ic[i]->Add(ilev,crack.c_old_mf);
211 }
212 }
213 else
214 {
215 crack.c_mf[ilev]->setVal(1.0);
216 crack.c_old_mf[ilev]->setVal(1.0);
217 }
218
219 if(material.is_ic) material.ic->Initialize(ilev,material.material_mf);
220 else material.material_mf[ilev]->setVal(1.0);
221 }
222
223 void TimeStepBegin(Set::Scalar /*time*/, int /*iter*/) override
224 {
225 Util::Message(INFO,crack.driving_force_norm," ",crack.driving_force_reference," ",crack.driving_force_tolerance_rel);
226 if (crack.driving_force_norm / crack.driving_force_reference < crack.driving_force_tolerance_rel)
227 elastic.do_solve_now = true;
228 if (crack.driving_force_norm < crack.driving_force_tolerance_abs)
229 elastic.do_solve_now = true;
230
231 if (!elastic.do_solve_now) return;
232 // if(iter%elastic.interval != 0) return;
233
234 Util::Message(INFO);
235 Util::Abort(INFO,"DegradeModulus is no longer implemented");
236 // for (int p = 0; p < number_of_materials; p++)
237 // material.brittlemodeltype[p].DegradeModulus(0.0);
238
239 for (int ilev = 0; ilev < nlevels; ++ilev)
240 {
241 //material.model_mf[ilev]->setVal(material.brittlemodeltype[0]);
242 Util::RealFillBoundary(*material.material_mf[ilev],geom[ilev]);
243 std::swap(elastic.energy_pristine_old_mf[ilev], elastic.energy_pristine_mf[ilev]);
244
245 for (MFIter mfi(*material.material_mf[ilev], false); mfi.isValid(); ++mfi)
246 {
247 amrex::Box bx = mfi.grownnodaltilebox();
248 amrex::Array4<brittle_fracture_model_type_test> const &model = material.model_mf[ilev]->array(mfi);
249 amrex::Array4<const Set::Scalar> const &mat = material.material_mf[ilev]->array(mfi);
250
251 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
252 brittle_fracture_model_type_test modelsum(0.0,0.0);
253 for (int p = 0; p < number_of_materials; p++)
254 {
255 modelsum += mat(i,j,k,p)*material.brittlemodeltype[p];
256 }
257 model(i,j,k) = /*(1./matsum) * */modelsum ;
258 });
259 }
260 Util::RealFillBoundary(*material.model_mf[ilev],geom[ilev]);
261 }
262
263 // Degrading the modulus
264 for (int ilev = 0; ilev < nlevels; ++ilev)
265 {
266 Util::RealFillBoundary(*crack.c_mf[ilev],geom[ilev]);
267 Util::RealFillBoundary(*material.material_mf[ilev],geom[ilev]);
268
269 for (amrex::MFIter mfi(*elastic.disp_mf[ilev],true); mfi.isValid(); ++mfi)
270 {
271 amrex::Box box = mfi.grownnodaltilebox();
272 amrex::Array4<const Set::Scalar> const& c = (*crack.c_mf[ilev]).array(mfi);
273 Util::Abort(INFO,"This is depricated and needs to be upgraded!");
274 //amrex::Array4<brittle_fracture_model_type_test> model = (material.model_mf)[ilev]->array(mfi);
275
276 amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k){
277 Set::Scalar _temp = 0.0;
278 // _temp = crack.cracktype->g_phi(c(i,j,k,0),0);
279 _temp = crack.cracktype.g_phi(c(i,j,k,0),0) + crack.scaleModulusMax;
280
281 if (std::isnan(_temp)) Util::Abort(INFO);
282 if(_temp < 0.0) _temp = 0.;
283 if(_temp > 1.0) _temp = 1.0;
284
285 // _temp = std::max(crack.scaleModulusMax, _temp);
286 // _temp = std::max(crack.scaleModulusMax, crack.scaleModulusMax + _temp * (1. - crack.scaleModulusMax));
287
288 Util::Abort(INFO,"DegradeModulus is no longer available");
289 //model(i,j,k).DegradeModulus(1-_temp);
290 });
291 }
292 Util::RealFillBoundary(*material.model_mf[ilev],geom[ilev]);
293 }
294
295 //
296 // Update the body force to be consistent with grid size
297 //
298 for (int ilev = 0; ilev < nlevels; ++ilev)
299 {
300 const Real* DX = geom[ilev].CellSize();
301 Set::Scalar volume = AMREX_D_TERM(DX[0],*DX[1],*DX[2]);
302
303 AMREX_D_TERM(elastic.rhs_mf[ilev]->setVal(loading.body_force(0)*volume,0,1);,
304 elastic.rhs_mf[ilev]->setVal(loading.body_force(1)*volume,1,1);,
305 elastic.rhs_mf[ilev]->setVal(loading.body_force(2)*volume,2,1););
306 }
307
308 elastic.brittlebc.Init(elastic.rhs_mf,geom);
309
310 //
311 // Prepare the operator for solve
312 //
313 Operator::Elastic<brittle_fracture_model_type_test::sym> op_b;
314 LPInfo info;
315 for (int ilev = 0; ilev < nlevels; ilev++) if (elastic.disp_mf[ilev]->contains_nan()) Util::Warning(INFO);
316
317 op_b.define(geom, grids, dmap, info);
318 op_b.SetBC(&elastic.brittlebc);
319 op_b.SetAverageDownCoeffs(true); // <<< Important! Do not change.
320 op_b.SetOmega(elastic.omega);
321 //
322 // Actually do the elastic solve
323 //
324 //for (int ilev = 0; ilev < nlevels; ilev++)
325 //{
326 // elastic.disp_mf[ilev]->setVal(0.0);
327 //}
328
329 {
330 IO::ParmParse pp("elastic");
331 Solver::Nonlocal::Newton<brittle_fracture_model_type_test> solver(op_b);
332 //Solver::Nonlocal::Linear<brittle_fracture_model_type_test> solver(op_b);
333 pp_queryclass("solver",solver);
334 solver.solve(elastic.disp_mf, elastic.rhs_mf, material.model_mf,1E-8,1E-8);
335 solver.compResidual(elastic.residual_mf,elastic.disp_mf,elastic.rhs_mf,material.model_mf);
336 }
337
338 for (int ilev = 0; ilev < nlevels; ilev++)
339 {
340 op_b.Strain(ilev,*elastic.strain_mf[ilev],*elastic.disp_mf[ilev]);
341 op_b.Stress(ilev,*elastic.stress_mf[ilev],*elastic.disp_mf[ilev]);
342 op_b.Energy(ilev,*elastic.energy_mf[ilev],*elastic.disp_mf[ilev]);
343
344 Util::RealFillBoundary(*elastic.strain_mf[ilev],geom[ilev]);
345 Util::RealFillBoundary(*elastic.energy_pristine_old_mf[ilev],geom[ilev]);
346 Util::RealFillBoundary(*material.material_mf[ilev],geom[ilev]);
347 elastic.energy_pristine_mf[ilev]->setVal(0.0);
348
349 for (amrex::MFIter mfi(*elastic.strain_mf[ilev],true); mfi.isValid(); ++mfi)
350 {
351 const amrex::Box& box = mfi.grownnodaltilebox();
352 amrex::Array4<const Set::Scalar> const& strain = (*elastic.strain_mf[ilev]).array(mfi);
353 amrex::Array4<const Set::Scalar> const& mat = (*material.material_mf[ilev]).array(mfi);
354 amrex::Array4<Set::Scalar> const& energy = (*elastic.energy_pristine_mf[ilev]).array(mfi);
355 amrex::Array4<Set::Scalar> const& energy_old = (*elastic.energy_pristine_old_mf[ilev]).array(mfi);
356
357 amrex::ParallelFor (box,[=] AMREX_GPU_DEVICE(int i, int j, int k){
358
359 Set::Matrix eps = Numeric::FieldToMatrix(strain,i,j,k);
360 Eigen::SelfAdjointEigenSolver<Set::Matrix> eigensolver(eps);
361 Set::Vector eValues = eigensolver.eigenvalues();
362 Set::Matrix eVectors = eigensolver.eigenvectors();
363
364 Set::Matrix epsp = Set::Matrix::Zero();
365 Set::Matrix epsn = Set::Matrix::Zero();
366
367 for (int n = 0; n < AMREX_SPACEDIM; n++)
368 {
369 if(eValues(n) > 0.0) epsp += eValues(n)*(eVectors.col(n)*eVectors.col(n).transpose());
370 else epsn += eValues(n)*(eVectors.col(n)*eVectors.col(n).transpose());
371 }
372 // brittle_fracture_model_type_test _tmpmodel(0.0,0.0);
373 // for (int p = 0; p < number_of_materials; p++)
374 // _tmpmodel += mat(i,j,k,p)*material.brittlemodeltype[p];
375
376
377 for (int p = 0; p < number_of_materials; p++)
378 energy(i,j,k) += mat(i,j,k,p)*material.brittlemodeltype[p].W(epsp);
379
380 if (energy(i,j,k) < energy_old(i,j,k)) energy(i,j,k) = energy_old(i,j,k);
381 });
382 }
383 Util::RealFillBoundary(*elastic.energy_pristine_mf[ilev],geom[ilev]);
384 }
385
386 integrate_variables_before_advance = false;
387 integrate_variables_after_advance = true;
388 }
389
390 void Advance(int lev, Set::Scalar /*time*/, Set::Scalar dt) override
391 {
392 std::swap(crack.c_old_mf[lev], crack.c_mf[lev]);
393
394 //Util::RealFillBoundary(*crack.c_old_mf[lev],geom[lev]);
395 //Util::RealFillBoundary(*crack.c_mf[lev],geom[lev]);
396 //Util::RealFillBoundary(*material.material_mf[lev],geom[lev]);
397 //Util::RealFillBoundary(*elastic.energy_pristine_mf[lev],geom[lev]);
398
399 const Set::Scalar* DX = geom[lev].CellSize();
400 amrex::Box domain(geom[lev].Domain());
401 domain.convert(amrex::IntVect::TheNodeVector());
402 const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
403
404 for ( amrex::MFIter mfi(*crack.c_mf[lev],true); mfi.isValid(); ++mfi )
405 {
406 const amrex::Box& bx = mfi.nodaltilebox();
407 amrex::Array4<const Set::Scalar> const& c_old = (*crack.c_old_mf[lev]).array(mfi);
408 amrex::Array4<Set::Scalar> const& df = (*crack.driving_force_mf[lev]).array(mfi);
409 amrex::Array4<Set::Scalar> const& c = (*crack.c_mf[lev]).array(mfi);
410 amrex::Array4<const Set::Scalar> const& energy = (*elastic.energy_pristine_mf[lev]).array(mfi);
411 amrex::Array4<const Set::Scalar> const& mat = (*material.material_mf[lev]).array(mfi);
412
413 amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k){
414
415#if AMREX_SPACEDIM !=2
416 Util::Abort(INFO, "This doesn't work with 1D or 3D yet");
417#endif
418 if (i == lo.x && j == lo.y) c(i,j,k,0) = c(i+1,j+1,k,0);
419 else if (i == lo.x && j == hi.y) c(i,j,k,0) = c(i+1,j-1,k,0);
420 else if (i == hi.x && j == lo.y) c(i,j,k,0) = c(i-1,j+1,k,0);
421 else if (i == hi.x && j == hi.y) c(i,j,k,0) = c(i-1,j-1,k,0);
422 else if (i == lo.x) c(i,j,k,0) = c(i+1,j,k,0);
423 else if (j == lo.y) c(i,j,k,0) = c(i,j+1,k,0);
424 else if (i == hi.x) c(i,j,k,0) = c(i-1,j,k,0);
425 else if (j == hi.y) c(i,j,k,0) = c(i,j-1,k,0);
426 else
427 {
428 Set::Scalar rhs = 0.0;
429 //Set::Vector Dc = Numeric::Gradient(c_old, i, j, k, 0, DX);
430
431 Set::Scalar bilap =
432 Numeric::Stencil<Set::Scalar, 4, 0, 0>::D(c_old,i,j,k,0,DX) +
433 2.0 * Numeric::Stencil<Set::Scalar, 2, 2, 0>::D(c_old,i,j,k,0,DX) +
434 Numeric::Stencil<Set::Scalar, 0, 4, 0>::D(c_old,i,j,k,0,DX);
435
436 Set::Scalar en_cell = energy(i,j,k,0);
437
438 //if (std::isnan(en_cell)) Util::Abort(INFO, "Nans detected in en_cell. energy_box(i,j,k,0) = ", energy(i,j,k,0));
439 //if (std::isinf(c_old(i,j,k,0))) Util::Abort(INFO, "Infs detected in c_old");
440 //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));
441
442 df(i,j,k,0) = crack.cracktype.Dg_phi(c_old(i,j,k,0),0.0)*en_cell*elastic.df_mult;
443 rhs += crack.cracktype.Dg_phi(c_old(i,j,k,0),0.0)*en_cell*elastic.df_mult;
444
445 //if(rhs > 1.e10) Util::Abort(INFO, "Very large values of rhs at (",i,",",j,",",k,") = ", rhs);
446
447 //Set::Matrix DDc = Numeric::Hessian(c_old, i, j, k, 0, DX);
448 Set::Scalar laplacian = Numeric::Laplacian(c_old,i,j,k,0,DX);//DDc.trace();
449
450 Set::Scalar Gc = crack.cracktype.Gc(0.0);
451 Set::Scalar Zeta = crack.cracktype.Zeta(0.0);
452 Set::Scalar _temp_product = 1.0;
453 for (int m = 0; m < number_of_materials; m++)
454 _temp_product *= mat(i,j,k,m);
455
456 if (number_of_materials > 1) Gc *= (1.0 - _temp_product*(1.-crack.mult_df_Gc)/0.25);
457
458 // Util::Message(INFO, "Gc = ", Gc);
459
460 df(i,j,k,1) = Gc*crack.cracktype.Dw_phi(c_old(i,j,k,0),0.0)/(4.0*Zeta);
461 rhs += Gc*crack.cracktype.Dw_phi(c_old(i,j,k,0),0.)/(4.0*Zeta);
462
463 df(i,j,k,2) = 2.0*Zeta*Gc*laplacian*crack.mult_df_lap;
464 rhs -= 2.0*Zeta*Gc*laplacian*crack.mult_df_lap;
465
466 df(i,j,k,3) = crack.beta*bilap;
467 rhs += crack.beta * bilap;
468
469 //if (!material.homogeneous)
470 //rhs *= crack.cracktype->g_phi(_temp,0.0);
471
472 //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));
473 df(i,j,k,4) = std::max(0.,rhs - crack.cracktype.DrivingForceThreshold(c_old(i,j,k,0)));
474 c(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));
475
476 if (c (i,j,k,0) < 0.0) c(i,j,k,0) = 0.0;
477 if (c (i,j,k,0) > 1.0) c(i,j,k,0) = 1.0;
478 }
479 });
480 }
481 Util::RealFillBoundary(*crack.c_mf[lev],geom[lev]);
482 }
483
484 void TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, amrex::Real /*time*/, int /*ngrow*/) override
485 {
486 const amrex::Real *DX = geom[lev].CellSize();
487 const Set::Vector dx(DX);
488 const Set::Scalar dxnorm = dx.lpNorm<2>();
489 amrex::Box domain(geom[lev].Domain());
490 domain.convert(amrex::IntVect::TheNodeVector());
491
492 for (amrex::MFIter mfi(*crack.c_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
493 {
494 const amrex::Box bx = mfi.tilebox() & domain.grow(-1);
495 amrex::Array4<char> const &tags = a_tags.array(mfi);
496 amrex::Array4<const Set::Scalar> const &c = (*crack.c_mf[lev]).array(mfi);
497 amrex::Array4<const Set::Scalar> const &mat = (*material.material_mf[lev]).array(mfi);
498
499 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
500 Set::Vector grad = Numeric::Gradient(c, i, j, k, 0, DX);
501 Set::Vector grad2 = Numeric::Gradient(mat, i, j, k, 0, DX);
502 if (dxnorm * grad.lpNorm<2>() >= crack.refinement_threshold || dxnorm * grad2.lpNorm<2>() >= material.refinement_threshold)
503 tags(i, j, k) = amrex::TagBox::SET;
504 });
505 }
506 }
507
508 void Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,const amrex::MFIter &mfi, const amrex::Box &box) override
509 {
510 const amrex::Real* DX = geom[amrlev].CellSize();
511 const Set::Scalar DV = AMREX_D_TERM(DX[0],*DX[1],*DX[2]);
512
513 amrex::Array4<const Set::Scalar> const &df = (*crack.driving_force_mf[amrlev]).array(mfi);
514 amrex::Array4<const Set::Scalar> const &c_new = (*crack.c_mf[amrlev]).array(mfi);
515 amrex::Array4<const Set::Scalar> const &energy = (*elastic.energy_mf[amrlev]).array(mfi);
516
517 amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
518 {
519 crack.driving_force_norm += df(i,j,k,4) * DV;
520 crack.int_crack += c_new(i,j,k) * DV;
521 elastic.int_energy += energy(i,j,k) * DV;
522 });
523 }
524
525 void TimeStepComplete(amrex::Real /*time*/,int /*iter*/) override
526 {
527 if (elastic.do_solve_now)
528 crack.driving_force_reference = crack.driving_force_norm;
529 //crack.driving_force_reference = crack.driving_force_norminf;
530
531 elastic.do_solve_now = false;
532 }
533
534private:
535
536 int number_of_ghost_nodes = 2; ///< Number of ghost nodes
537 int number_of_materials = 1;
538 int nlevels;
539
540 struct{
541 Set::Field<Set::Scalar> disp_mf; ///< displacement field
542 Set::Field<Set::Scalar> strain_mf; ///< total strain field (gradient of displacement)
543 Set::Field<Set::Scalar> stress_mf; ///< stress field
544 Set::Field<Set::Scalar> rhs_mf; ///< rhs fab for elastic solution
545 Set::Field<Set::Scalar> residual_mf; ///< residual field for solver
546 Set::Field<Set::Scalar> energy_mf; ///< total elastic energy
547 Set::Field<Set::Scalar> energy_pristine_mf; ///< energy of the prisitne material as if no crack is present
548 Set::Field<Set::Scalar> energy_pristine_old_mf; ///< energy of the pristine material for previous time step.
549 Set::Scalar int_energy = 0.0; ///< integrated energy over the entire domain.
550
551 BC::Operator::Elastic::Constant brittlebc; ///< elastic BC if using brittle fracture
552 Set::Scalar df_mult = 1.0; ///< mulitplier for elastic driving force.
553 bool do_solve_now = false;
554 int interval = 0;
555 Set::Scalar omega = 2./3.;
556 } elastic;
557
558
559
560 struct{
561 Set::Field<Set::Scalar> c_mf; ///< crack field at current time step
562 Set::Field<Set::Scalar> c_old_mf; ///< crack field at previous time step
563 Set::Field<Set::Scalar> driving_force_mf; ///< crack driving forces.
564 Set::Scalar int_crack = 0.0; ///< integrated crack field over the entire domain.
565 Set::Scalar driving_force_reference = 1.0;
566 Set::Scalar driving_force_norm = 0.0;
567 Set::Scalar driving_force_norminf = 0.0;
568 Set::Scalar driving_force_tolerance_rel = 1E-4;
569 Set::Scalar driving_force_tolerance_abs = 0.0;
570
571 Model::Interface::Crack::Constant cracktype; ///< type of crack. See Crack/Constant or Crack/Sin
572 amrex::Vector<std::string> ic_type; ///< crack IC type. See IC/Notch and IC/Ellipsoid
573 amrex::Vector<IC::IC<Set::Scalar>*> ic; ///< crack IC. See IC/Notch and IC/Ellipsoid
574 bool is_ic = false;
575
576 Set::Scalar scaleModulusMax = 0.02; ///< material modulus ratio inside crack (default = 0.02).
577 Set::Scalar refinement_threshold = 0.001;
578
579 Set::Scalar mult_df_Gc = 1.0; ///< Multiplier for Gc/zeta term
580 Set::Scalar mult_df_lap = 1.0; ///< Multiplier for the laplacian term
581 Set::Scalar beta = 0.0;
582 } crack;
583
584 struct{
585 amrex::Vector<brittle_fracture_model_type_test> brittlemodeltype;
586 Set::Field<Set::Scalar> material_mf;
587 Set::Field<brittle_fracture_model_type_test> model_mf;
588 std::string input_material = "isotropic";
589 Set::Scalar refinement_threshold = 0.1;
590
591 IC::IC<Set::Scalar> *ic;
592 std::string ic_type;
593 bool is_ic = false;
594 } material;
595
596 struct{
597 Set::Vector body_force = Set::Vector::Zero();
598 Set::Scalar val = 0.;
599 } loading;
600
601};
602}
603#endif
pp_queryarr
#define pp_queryarr(...)
Definition ParmParse.H:103
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::contains
bool contains(std::string name)
Definition ParmParse.H:154
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::Fracture::material
struct Integrator::Fracture::@9 material
Integrator::Fracture::scaleModulusMax
Set::Scalar scaleModulusMax
material modulus ratio inside crack (default = 0.02).
Definition Fracture.H:576
Integrator::Fracture::stress_mf
Set::Field< Set::Scalar > stress_mf
stress field
Definition Fracture.H:543
Integrator::Fracture::ic
amrex::Vector< IC::IC< Set::Scalar > * > ic
crack IC. See IC/Notch and IC/Ellipsoid
Definition Fracture.H:573
Integrator::Fracture::crack
struct Integrator::Fracture::@8 crack
Integrator::Fracture::ic_type
amrex::Vector< std::string > ic_type
crack IC type. See IC/Notch and IC/Ellipsoid
Definition Fracture.H:572
Integrator::Fracture::ic
IC::IC< Set::Scalar > * ic
Definition Fracture.H:591
Integrator::Fracture::elastic
struct Integrator::Fracture::@7 elastic
Integrator::Fracture::rhs_mf
Set::Field< Set::Scalar > rhs_mf
rhs fab for elastic solution
Definition Fracture.H:544
Integrator::Fracture::driving_force_tolerance_abs
Set::Scalar driving_force_tolerance_abs
Definition Fracture.H:569
Integrator::Fracture::int_crack
Set::Scalar int_crack
integrated crack field over the entire domain.
Definition Fracture.H:564
Integrator::Fracture::Initialize
void Initialize(int ilev) override
Definition Fracture.H:191
Integrator::Fracture::Integrate
void Integrate(int amrlev, Set::Scalar, int, const amrex::MFIter &mfi, const amrex::Box &box) override
Definition Fracture.H:508
Integrator::Fracture::c_old_mf
Set::Field< Set::Scalar > c_old_mf
crack field at previous time step
Definition Fracture.H:562
Integrator::Fracture::loading
struct Integrator::Fracture::@10 loading
Integrator::Fracture::energy_pristine_old_mf
Set::Field< Set::Scalar > energy_pristine_old_mf
energy of the pristine material for previous time step.
Definition Fracture.H:548
Integrator::Fracture::Advance
void Advance(int lev, Set::Scalar, Set::Scalar dt) override
Definition Fracture.H:390
Integrator::Fracture::TimeStepBegin
void TimeStepBegin(Set::Scalar, int) override
Definition Fracture.H:223
Integrator::Fracture::energy_pristine_mf
Set::Field< Set::Scalar > energy_pristine_mf
energy of the prisitne material as if no crack is present
Definition Fracture.H:547
Integrator::Fracture::disp_mf
Set::Field< Set::Scalar > disp_mf
displacement field
Definition Fracture.H:541
Integrator::Fracture::brittlebc
BC::Operator::Elastic::Constant brittlebc
elastic BC if using brittle fracture
Definition Fracture.H:551
Integrator::Fracture::omega
Set::Scalar omega
Definition Fracture.H:555
Integrator::Fracture::TagCellsForRefinement
void TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, amrex::Real, int) override
Definition Fracture.H:484
Integrator::Fracture::body_force
Set::Vector body_force
Definition Fracture.H:597
Integrator::Fracture::number_of_ghost_nodes
int number_of_ghost_nodes
Number of ghost nodes.
Definition Fracture.H:536
Integrator::Fracture::cracktype
Model::Interface::Crack::Constant cracktype
type of crack. See Crack/Constant or Crack/Sin
Definition Fracture.H:571
Integrator::Fracture::mult_df_Gc
Set::Scalar mult_df_Gc
Multiplier for Gc/zeta term.
Definition Fracture.H:579
Integrator::Fracture::residual_mf
Set::Field< Set::Scalar > residual_mf
residual field for solver
Definition Fracture.H:545
Integrator::Fracture::name
static constexpr const char * name
Definition Fracture.H:48
Integrator::Fracture::input_material
std::string input_material
Definition Fracture.H:588
Integrator::Fracture::driving_force_norminf
Set::Scalar driving_force_norminf
Definition Fracture.H:567
Integrator::Fracture::mult_df_lap
Set::Scalar mult_df_lap
Multiplier for the laplacian term.
Definition Fracture.H:580
Integrator::Fracture::driving_force_tolerance_rel
Set::Scalar driving_force_tolerance_rel
Definition Fracture.H:568
Integrator::Fracture::ic_type
std::string ic_type
Definition Fracture.H:592
Integrator::Fracture::energy_mf
Set::Field< Set::Scalar > energy_mf
total elastic energy
Definition Fracture.H:546
Integrator::Fracture::model_mf
Set::Field< brittle_fracture_model_type_test > model_mf
Definition Fracture.H:587
Integrator::Fracture::TimeStepComplete
void TimeStepComplete(amrex::Real, int) override
Definition Fracture.H:525
Integrator::Fracture::df_mult
Set::Scalar df_mult
mulitplier for elastic driving force.
Definition Fracture.H:552
Integrator::Fracture::refinement_threshold
Set::Scalar refinement_threshold
Definition Fracture.H:577
Integrator::Fracture::brittlemodeltype
amrex::Vector< brittle_fracture_model_type_test > brittlemodeltype
Definition Fracture.H:585
Integrator::Fracture::material_mf
Set::Field< Set::Scalar > material_mf
Definition Fracture.H:586
Integrator::Fracture::driving_force_reference
Set::Scalar driving_force_reference
Definition Fracture.H:565
Integrator::Fracture::driving_force_mf
Set::Field< Set::Scalar > driving_force_mf
crack driving forces.
Definition Fracture.H:563
Integrator::Fracture::val
Set::Scalar val
Definition Fracture.H:598
Integrator::Fracture::int_energy
Set::Scalar int_energy
integrated energy over the entire domain.
Definition Fracture.H:549
Integrator::Fracture::strain_mf
Set::Field< Set::Scalar > strain_mf
total strain field (gradient of displacement)
Definition Fracture.H:542
Integrator::Fracture::c_mf
Set::Field< Set::Scalar > c_mf
crack field at current time step
Definition Fracture.H:561
Integrator::Fracture::driving_force_norm
Set::Scalar driving_force_norm
Definition Fracture.H:566
Integrator::Fracture::beta
Set::Scalar beta
Definition Fracture.H:581
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
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::SetAverageDownCoeffs
virtual void SetAverageDownCoeffs(bool a_average_down_coeffs) override
Definition Elastic.H:110
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
Operator::Operator< Grid::Node >::SetOmega
void SetOmega(Set::Scalar a_omega)
Definition Operator.H:57
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::brittle_fracture_model_type_test
Model::Solid::Linear::Isotropic brittle_fracture_model_type_test
Definition Fracture.H:41
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::Laplacian
AMREX_FORCE_INLINE Set::Scalar Laplacian(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:546
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