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