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 
39 namespace Integrator
40 {
41 //using brittle_fracture_model_type_test = Model::Solid::Linear::Isotropic;
42 using suture_fracture_model_type = Model::Solid::Linear::IsotropicDegradable;
43 
44 class SutureCrack : public Integrator
45 {
46 
47 public:
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 
148 protected:
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 
465 private:
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
Integrator::SutureCrack::verbose
int verbose
Definition: SutureCrack.H:527
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
Integrator::SutureCrack::max_iter
int max_iter
Definition: SutureCrack.H:523
Integrator::SutureCrack
Definition: SutureCrack.H:44
Integrator::SutureCrack::linop_maxorder
int linop_maxorder
Definition: SutureCrack.H:536
Integrator::SutureCrack::Initialize
void Initialize(int ilev) override
Definition: SutureCrack.H:149
Integrator::SutureCrack::agglomeration
bool agglomeration
Definition: SutureCrack.H:539
Util::RealFillBoundary
AMREX_FORCE_INLINE void RealFillBoundary(amrex::FabArray< amrex::BaseFab< T >> &a_mf, const amrex::Geometry &)
Definition: Util.H:307
IC::IC
Pure abstract IC object from which all other IC objects inherit.
Definition: IC.H:12
Integrator::SutureCrack::material
struct Integrator::SutureCrack::@24 material
Operator::Elastic
Definition: Elastic.H:22
Integrator::Integrator::box
std::vector< amrex::Box > box
Definition: Integrator.H:367
Integrator::SutureCrack::mult_df_lap
Set::Scalar mult_df_lap
Multiplier for the laplacian term.
Definition: SutureCrack.H:504
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::do_solve_now
bool do_solve_now
Definition: SutureCrack.H:483
Integrator::SutureCrack::field_old
Set::Field< Set::Scalar > field_old
crack field at previous time step
Definition: SutureCrack.H:488
Integrator::SutureCrack::refinement_threshold
Set::Scalar refinement_threshold
Definition: SutureCrack.H:501
Integrator::SutureCrack::Advance
void Advance(int lev, Set::Scalar, Set::Scalar dt) override
Definition: SutureCrack.H:330
Integrator::SutureCrack::val
Set::Scalar val
Definition: SutureCrack.H:547
Solver::Nonlocal::Linear::setVerbose
void setVerbose(const int a_verbose)
Definition: Linear.H:184
Integrator::SutureCrack::brittlebc
BC::Operator::Elastic::Constant brittlebc
elastic BC if using brittle fracture
Definition: SutureCrack.H:481
Integrator.H
Notch.H
Set::Field< Set::Scalar >
Definition: Set.H:236
Solver::Nonlocal::Linear::setMaxIter
void setMaxIter(const int a_max_iter)
Definition: Linear.H:180
Integrator::SutureCrack::rhs
Set::Field< Set::Scalar > rhs
rhs fab for elastic solution
Definition: SutureCrack.H:475
Integrator::suture_fracture_model_type
Model::Solid::Linear::IsotropicDegradable suture_fracture_model_type
Definition: SutureCrack.H:42
Integrator::Integrator::integrate_variables_after_advance
bool integrate_variables_after_advance
Definition: Integrator.H:314
Laminate.H
Solver::Nonlocal::Newton
Definition: Newton.H:15
Integrator::SutureCrack::modulus_field
Set::Field< Set::Scalar > modulus_field
Definition: SutureCrack.H:509
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition: Base.H:20
Constant.H
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:655
pp_query
#define pp_query(...)
Definition: ParmParse.H:104
Integrator::SutureCrack::ic
IC::IC * ic
crack IC. See IC/Notch and IC/Ellipsoid
Definition: SutureCrack.H:497
Solver::Nonlocal::Linear::setPreSmooth
void setPreSmooth(const int a_pre_smooth)
Definition: Linear.H:185
Integrator::SutureCrack::residual
Set::Field< Set::Scalar > residual
residual field for solver
Definition: SutureCrack.H:476
Elastic.H
Integrator::SutureCrack::tol_rel
Set::Scalar tol_rel
Definition: SutureCrack.H:529
Integrator::SutureCrack::driving_force
Set::Field< Set::Scalar > driving_force
crack driving forces.
Definition: SutureCrack.H:489
Integrator::SutureCrack::cg_tol_abs
Set::Scalar cg_tol_abs
Definition: SutureCrack.H:532
Set::Scalar
amrex::Real Scalar
Definition: Base.H:19
Integrator::SutureCrack::driving_force_tolerance_rel
Set::Scalar driving_force_tolerance_rel
Definition: SutureCrack.H:492
Integrator::SutureCrack::cgverbose
int cgverbose
Definition: SutureCrack.H:528
Integrator::SutureCrack::tend
Set::Scalar tend
Definition: SutureCrack.H:534
IC::Notch
Definition: Notch.H:17
pp_queryclass
#define pp_queryclass(...)
Definition: ParmParse.H:105
Solid.H
Integrator::SutureCrack::strain
Set::Field< Set::Scalar > strain
total strain field (gradient of displacement)
Definition: SutureCrack.H:473
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition: Base.H:23
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
Integrator::SutureCrack::loading
struct Integrator::SutureCrack::@26 loading
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
Sin.H
Crack.H
Integrator::SutureCrack::nlevels
int nlevels
Definition: SutureCrack.H:469
Integrator::SutureCrack::Integrate
void Integrate(int amrlev, Set::Scalar, int, const amrex::MFIter &mfi, const amrex::Box &box)
Definition: SutureCrack.H:444
Integrator::SutureCrack::interval
int interval
Definition: SutureCrack.H:521
Solver::Nonlocal::Linear::setPostSmooth
void setPostSmooth(const int a_post_smooth)
Definition: Linear.H:186
Integrator::Integrator::RegisterGeneralFab
void RegisterGeneralFab(Set::Field< T > &new_fab, int ncomp, int nghost, bool evolving=true)
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::SutureCrack::bottom_tol
Set::Scalar bottom_tol
Definition: SutureCrack.H:520
BC::Operator::Elastic::Constant
Definition: Constant.H:31
Integrator::SutureCrack::sol
struct Integrator::SutureCrack::@25 sol
Integrator::SutureCrack::number_of_materials
int number_of_materials
Definition: SutureCrack.H:468
Integrator::SutureCrack::bottom_max_iter
int bottom_max_iter
Definition: SutureCrack.H:525
Constant.H
Integrator::SutureCrack::driving_force_tolerance_abs
Set::Scalar driving_force_tolerance_abs
Definition: SutureCrack.H:493
Integrator::SutureCrack::type
std::string type
Definition: SutureCrack.H:522
Integrator::SutureCrack::stress
Set::Field< Set::Scalar > stress
stress field
Definition: SutureCrack.H:474
Integrator::SutureCrack::is_ic
bool is_ic
Definition: SutureCrack.H:498
Model::Interface::Crack::Constant
Definition: Constant.H:18
Integrator::SutureCrack::disp
Set::Field< Set::Scalar > disp
displacement field
Definition: SutureCrack.H:472
Integrator::SutureCrack::input_material
std::string input_material
Definition: SutureCrack.H:511
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:165
Integrator
Collection of numerical integrator objects.
Definition: AllenCahn.H:39
Integrator::SutureCrack::SutureCrack
SutureCrack()
Definition: SutureCrack.H:48
Integrator::SutureCrack::mult_df_Gc
Set::Scalar mult_df_Gc
Multiplier for Gc/zeta term.
Definition: SutureCrack.H:503
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
Integrator::SutureCrack::driving_force_reference
Set::Scalar driving_force_reference
Definition: SutureCrack.H:490
Integrator::SutureCrack::number_of_ghost_nodes
int number_of_ghost_nodes
Number of ghost nodes.
Definition: SutureCrack.H:467
Newton.H
Integrator::SutureCrack::brittlemodel
Set::Field< suture_fracture_model_type > brittlemodel
Definition: SutureCrack.H:510
Integrator::SutureCrack::scaleModulusMax
Set::Scalar scaleModulusMax
material modulus ratio inside crack (default = 0.02).
Definition: SutureCrack.H:500
Integrator::SutureCrack::crack
struct Integrator::SutureCrack::@23 crack
Integrator::SutureCrack::bottom_solver
std::string bottom_solver
Definition: SutureCrack.H:535
Operator::Elastic::Strain
void Strain(int amrlev, amrex::MultiFab &epsfab, const amrex::MultiFab &ufab, bool voigt=false) const
Definition: Elastic.cpp:403
Operator::Elastic::SetBC
void SetBC(::BC::Operator::Elastic::Elastic *a_bc)
Definition: Elastic.H:61
Operator::Elastic::Energy
void Energy(int amrlev, amrex::MultiFab &energy, const amrex::MultiFab &u, bool a_homogeneous=false)
Definition: Elastic.cpp:529
Integrator::SutureCrack::TagCellsForRefinement
void TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, amrex::Real, int) override
Definition: SutureCrack.H:420
Integrator::Integrator::RegisterNodalFab
void RegisterNodalFab(Set::Field< Set::Scalar > &new_fab, int ncomp, int nghost, std::string name, bool writeout)
Definition: Integrator.cpp:310
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:909
Integrator::SutureCrack::consolidation
bool consolidation
Definition: SutureCrack.H:540
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:1027
Integrator::SutureCrack::body_force
Set::Vector body_force
Definition: SutureCrack.H:546
Util::Warning
void Warning(std::string file, std::string func, int line, Args const &... args)
Definition: Util.H:171
Integrator::SutureCrack::max_fixed_iter
int max_fixed_iter
Definition: SutureCrack.H:526
BC.H
IO::ParmParse
Definition: ParmParse.H:110
Integrator::SutureCrack::TimeStepComplete
void TimeStepComplete(amrex::Real, int)
Definition: SutureCrack.H:457
Integrator::SutureCrack::tstart
Set::Scalar tstart
Definition: SutureCrack.H:533
Integrator::SutureCrack::TimeStepBegin
void TimeStepBegin(Set::Scalar, int iter) override
Definition: SutureCrack.H:178
Integrator::SutureCrack::energy
Set::Field< Set::Scalar > energy
total elastic energy
Definition: SutureCrack.H:477
Integrator::SutureCrack::field
Set::Field< Set::Scalar > field
crack field at current time step
Definition: SutureCrack.H:487
Integrator::SutureCrack::cracktype
Model::Interface::Crack::Crack * cracktype
type of crack. See Crack/Constant or Crack/Sin
Definition: SutureCrack.H:495
Integrator::SutureCrack::elastic
struct Integrator::SutureCrack::@22 elastic
Operator::Elastic::Stress
void Stress(int amrlev, amrex::MultiFab &sigmafab, const amrex::MultiFab &ufab, bool voigt=false, bool a_homogeneous=false)
Definition: Elastic.cpp:464
Integrator::SutureCrack::use_fsmooth
bool use_fsmooth
Definition: SutureCrack.H:537
Integrator::SutureCrack::pre_smooth
int pre_smooth
Definition: SutureCrack.H:541
Integrator::SutureCrack::tol_abs
Set::Scalar tol_abs
Definition: SutureCrack.H:530
Solver::Nonlocal::Linear::setBottomMaxIter
void setBottomMaxIter(const int a_bottom_max_iter)
Definition: Linear.H:181
Solver::Nonlocal::Linear::setFixedIter
void setFixedIter(const int a_fixed_iter)
Definition: Linear.H:183
IO::ParmParse::queryclass
void queryclass(std::string name, T *value, std::string file="", std::string func="", int line=-1)
Definition: ParmParse.H:443
Integrator::Integrator::integrate_variables_before_advance
bool integrate_variables_before_advance
Definition: Integrator.H:313
Integrator::SutureCrack::brittlemodeltype
suture_fracture_model_type brittlemodeltype
Definition: SutureCrack.H:508
Integrator::SutureCrack::post_smooth
int post_smooth
Definition: SutureCrack.H:542
IC::Laminate
Definition: Laminate.H:16
Integrator::SutureCrack::max_fmg_iter
int max_fmg_iter
Definition: SutureCrack.H:524
Isotropic.H
Solver::Nonlocal::Linear::setMaxFmgIter
void setMaxFmgIter(const int a_max_fmg_iter)
Definition: Linear.H:182
INFO
#define INFO
Definition: Util.H:20
Integrator::Integrator::dt
amrex::Vector< amrex::Real > dt
Timesteps for each level of refinement.
Definition: Integrator.H:303
IC.H
Set::Field< suture_fracture_model_type >
Integrator::SutureCrack::cg_tol_rel
Set::Scalar cg_tol_rel
Definition: SutureCrack.H:531
IO::ParmParse::queryarr
int queryarr(std::string name, std::vector< T > &value, std::string, std::string, int)
Definition: ParmParse.H:333
Integrator::SutureCrack::max_coarsening_level
int max_coarsening_level
Definition: SutureCrack.H:538
Integrator::Integrator::RegisterIntegratedVariable
void RegisterIntegratedVariable(Set::Scalar *integrated_variable, std::string name, bool extensive=true)
Definition: Integrator.cpp:320
Integrator::SutureCrack::driving_force_norm
Set::Scalar driving_force_norm
Definition: SutureCrack.H:491
Linear.H
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition: Util.H:133
Model::Interface::Crack::Crack
Definition: Crack.H:16
Generated on Fri Aug 16 2024 16:30:29 for Alamo by   doxygen 1.8.17