Line data Source code
1 : //
2 : // Create a simple notch in an otherwise uniformly filled region.
3 : // (This was created for, and mostly used for, Mode II fracture tests.)
4 : //
5 : // This is an old IC that should be replaced by IC::Expression
6 : //
7 :
8 : #ifndef IC_NOTCH_H_
9 : #define IC_NOTCH_H_
10 :
11 : #include "IC/IC.H"
12 : #include "IO/ParmParse.H"
13 : #include "Set/Set.H"
14 :
15 : // Note: right now this is meant for 2D. We need to rethink this implementation for 3D.
16 : namespace IC
17 : {
18 : class Notch : public IC<Set::Scalar>
19 : {
20 : public:
21 : static constexpr const char *name = "notch";
22 :
23 : public:
24 : enum Mollifier
25 : {
26 : Dirac,
27 : Erf,
28 : Gaussian,
29 : Cosine,
30 : Linear
31 : };
32 :
33 : Notch(amrex::Vector<amrex::Geometry> &_geom) : IC(_geom)
34 : {
35 : nCenter.resize(1);
36 : nCenter[0] = Set::Vector::Zero();
37 : nOrientation.resize(1);
38 : nOrientation[0] = Set::Vector::Random();
39 : nThickness.resize(1);
40 : nThickness[0] = 0.01;
41 : nLength.resize(1);
42 : nLength[0] = 0.1;
43 : moll = Mollifier::Dirac;
44 : eps = 1e-5;
45 : }
46 0 : Notch(amrex::Vector<amrex::Geometry> &_geom, IO::ParmParse &pp, std::string name) : IC<Set::Scalar>(_geom)
47 : {
48 0 : pp_queryclass(name, *this);
49 0 : }
50 :
51 : void
52 0 : Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
53 : {
54 0 : Set::Vector DX(geom[lev].CellSize());
55 0 : Set::Scalar pi = std::atan(1.0) * 4;
56 :
57 0 : for (amrex::MFIter mfi(*a_field[lev], true); mfi.isValid(); ++mfi)
58 : {
59 : // amrex::Box bx = mfi.grownnodaltilebox();
60 0 : amrex::Box bx = mfi.tilebox();
61 0 : bx.grow(a_field[lev]->nGrow());
62 0 : amrex::Array4<Set::Scalar> const &field = a_field[lev]->array(mfi);
63 0 : amrex::IndexType type = a_field[lev]->ixType();
64 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
65 0 : Set::Vector x;
66 : // NODE
67 0 : if (type == amrex::IndexType::TheNodeType())
68 : {
69 0 : AMREX_D_TERM( x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i)) * geom[lev].CellSize()[0];,
70 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j)) * geom[lev].CellSize()[1];,
71 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k)) * geom[lev].CellSize()[2];);
72 : }
73 0 : else if (type == amrex::IndexType::TheCellType())
74 : {
75 0 : AMREX_D_TERM( x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
76 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
77 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
78 : }
79 :
80 0 : Set::Scalar min_value = field(i, j, k);
81 :
82 0 : for (int m = 0; m < nCenter.size(); m++)
83 : {
84 0 : Set::Scalar value = 1.0;
85 0 : if (std::abs((x - nCenter[m]).transpose() * nOrientation[m]) <= 0.5 * nLength[m])
86 : {
87 0 : Set::Scalar t = std::abs((x - nCenter[m]).transpose() * nNormal[m]) - (nThickness[m] / 2.0);
88 0 : if (moll == Mollifier::Erf)
89 0 : value = 0.5 + 0.5 * std::erf(t / eps);
90 : // value = std::erf(t/eps);
91 0 : else if (moll == Mollifier::Gaussian)
92 : {
93 0 : value = (t >= 0.0) ? 1.0 - std::exp(-0.5 * (4.0 * t / eps) * (4.0 * t / eps)) : 0.0;
94 : }
95 0 : else if (moll == Mollifier::Linear)
96 : {
97 0 : value = (t > eps) ? 1.0 : t / eps;
98 : }
99 0 : else if (moll == Mollifier::Dirac)
100 : {
101 0 : if (t < 0.0)
102 0 : value = 0.0;
103 : else
104 0 : value = 1.0;
105 : }
106 0 : else if (moll == Mollifier::Cosine)
107 : {
108 0 : if (t < 0.0)
109 0 : value = 0.0;
110 0 : else if (t < eps)
111 0 : value = 0.5 - 0.5 * std::cos(pi * t / eps);
112 : else
113 0 : value = 1.0;
114 : }
115 : }
116 : else
117 : {
118 0 : Set::Vector nLeft = nCenter[m] - 0.5 * nLength[m] * nOrientation[m];
119 0 : Set::Vector nRight = nCenter[m] + 0.5 * nLength[m] * nOrientation[m];
120 0 : Set::Vector nLeft2 = nLeft, nRight2 = nRight;
121 0 : Set::Vector correction = nOrientation[m] * (std::sqrt(nRadius[m] * nRadius[m] - nThickness[m] * nThickness[m] / 4.0));
122 0 : if (nRadius[m] > 0.5 * nThickness[m])
123 : {
124 0 : nLeft2 = nLeft + correction;
125 0 : nRight2 = nRight - correction;
126 : }
127 :
128 0 : Set::Scalar cosR = (nRadius[m] <= 0.0) ? 1.0 : std::sqrt(nRadius[m] * nRadius[m] - nThickness[m] * nThickness[m] / 4.0) / nRadius[m];
129 0 : Set::Scalar cosLeft = (x - nLeft2).dot(nOrientation[m]) / ((x - nLeft2).lpNorm<2>());
130 0 : Set::Scalar sinLeft = (x - nLeft2).dot(nNormal[m]) / ((x - nLeft2).lpNorm<2>());
131 0 : Set::Scalar cosRight = (x - nRight2).dot(nOrientation[m]) / ((x - nRight2).lpNorm<2>());
132 0 : Set::Scalar sinRight = (x - nRight2).dot(nNormal[m]) / ((x - nRight2).lpNorm<2>());
133 :
134 0 : Set::Scalar distLeft = (x - nLeft).dot(nOrientation[m]);
135 0 : Set::Scalar distRight = (x - nRight).dot(nOrientation[m]);
136 :
137 0 : if (distLeft <= 0.)
138 : {
139 0 : Set::Scalar t = (x - nLeft2).lpNorm<2>() - nRadius[m];
140 0 : if (moll == Mollifier::Erf)
141 : {
142 0 : value = 0.5 + 0.5 * std::erf(((x - nLeft).lpNorm<2>() - nThickness[m] / 2.0) / eps);
143 : }
144 0 : else if (moll == Mollifier::Gaussian)
145 : {
146 0 : t = ((x - nLeft).lpNorm<2>() - nThickness[m] / 2.0);
147 0 : value = (t >= 0.0) ? 1.0 - std::exp(-0.5 * (4.0 * t / eps) * (4.0 * t / eps)) : 0.0;
148 : }
149 0 : else if (moll == Mollifier::Linear)
150 : {
151 0 : t = ((x - nLeft).lpNorm<2>() - nThickness[m] / 2.0);
152 0 : value = (t >= 0.0) ? ((t > eps) ? 1.0 : t / eps) : 0.0;
153 : }
154 0 : else if (moll == Mollifier::Dirac)
155 : {
156 0 : if (t < 0.0)
157 0 : value = 0.0;
158 : else
159 0 : value = 1.0;
160 : }
161 0 : else if (moll == Mollifier::Cosine)
162 : {
163 0 : Set::Vector nTopLeft = nLeft + 0.5 * nThickness[m] * nNormal[m];
164 0 : Set::Vector nBotLeft = nLeft - 0.5 * nThickness[m] * nNormal[m];
165 :
166 0 : if (nRadius[m] <= 0.0)
167 : {
168 0 : Set::Scalar t2 = -(x - nTopLeft).dot(nOrientation[m]);
169 :
170 0 : Set::Scalar angTop = (x - nTopLeft).dot(nNormal[m]);
171 0 : Set::Scalar angBot = (x - nBotLeft).dot(nNormal[m]);
172 :
173 0 : if (angTop <= 0.0 && angBot >= 0.0)
174 : {
175 0 : if (t2 < 0.0)
176 0 : value = 0.0;
177 0 : else if (t2 < eps)
178 0 : value = 0.5 - 0.5 * std::cos(pi * t2 / eps);
179 : else
180 0 : value = 1.0;
181 : }
182 0 : else if (angTop > 0.0 && angBot > 0.0)
183 : {
184 0 : Set::Scalar t3 = (x - nTopLeft).lpNorm<2>();
185 0 : if (t3 < 0.0)
186 0 : value = 0.0;
187 0 : else if (t3 < eps)
188 0 : value = 0.5 - 0.5 * std::cos(pi * t3 / eps);
189 : else
190 0 : value = 1.0;
191 0 : }
192 : else
193 : {
194 0 : Set::Scalar t3 = (x - nBotLeft).lpNorm<2>();
195 0 : if (t3 < 0.0)
196 0 : value = 0.0;
197 0 : else if (t3 < eps)
198 0 : value = 0.5 - 0.5 * std::cos(pi * t3 / eps);
199 : else
200 0 : value = 1.0;
201 : }
202 : }
203 :
204 : else
205 : {
206 0 : if (std::abs(cosLeft) > cosR)
207 : {
208 0 : if (t < 0.0)
209 0 : value = 0.0;
210 0 : else if (t < eps)
211 0 : value = 0.5 - 0.5 * std::cos(pi * t / eps);
212 : else
213 0 : value = 1.0;
214 : }
215 :
216 : else
217 : {
218 0 : if (sinLeft >= 0)
219 : {
220 0 : Set::Scalar t2 = (x - nTopLeft).lpNorm<2>();
221 0 : if (t2 < 0.0)
222 0 : value = 0.0;
223 0 : else if (t2 < eps)
224 0 : value = 0.5 - 0.5 * std::cos(pi * t2 / eps);
225 : else
226 0 : value = 1.0;
227 : }
228 : else
229 : {
230 0 : Set::Scalar t2 = (x - nBotLeft).lpNorm<2>();
231 0 : if (t2 < 0.0)
232 0 : value = 0.0;
233 0 : else if (t2 < eps)
234 0 : value = 0.5 - 0.5 * std::cos(pi * t2 / eps);
235 : else
236 0 : value = 1.0;
237 : }
238 : }
239 : }
240 : }
241 : }
242 0 : if (distRight >= 0.)
243 : {
244 0 : Set::Scalar t = (x - nRight2).lpNorm<2>() - nRadius[m];
245 0 : if (moll == Mollifier::Erf)
246 : {
247 0 : value = 0.5 + 0.5 * std::erf(((x - nRight).lpNorm<2>() - nThickness[m] / 2.0) / eps);
248 : }
249 0 : else if (moll == Mollifier::Gaussian)
250 : {
251 0 : t = ((x - nRight).lpNorm<2>() - nThickness[m] / 2.0);
252 0 : value = (t >= 0.0) ? 1.0 - std::exp(-0.5 * (4.0 * t / eps) * (4.0 * t / eps)) : 0.0;
253 : }
254 0 : else if (moll == Mollifier::Linear)
255 : {
256 0 : t = ((x - nRight).lpNorm<2>() - nThickness[m] / 2.0);
257 0 : value = (t >= 0.0) ? ((t > eps) ? 1.0 : t / eps) : 0.0;
258 : }
259 0 : else if (moll == Mollifier::Dirac)
260 : {
261 0 : if (t < 0.0)
262 0 : value = 0.0;
263 : else
264 0 : value = 1.0;
265 : }
266 0 : else if (moll == Mollifier::Cosine)
267 : {
268 0 : Set::Vector nTopRight = nRight + 0.5 * nThickness[m] * nNormal[m];
269 0 : Set::Vector nBotRight = nRight - 0.5 * nThickness[m] * nNormal[m];
270 :
271 0 : if (nRadius[m] <= 0.0)
272 : {
273 0 : Set::Scalar t2 = (x - nTopRight).dot(nOrientation[m]);
274 0 : Set::Scalar angTop = (x - nTopRight).dot(nNormal[m]);
275 0 : Set::Scalar angBot = (x - nBotRight).dot(nNormal[m]);
276 :
277 0 : if (angTop <= 0.0 && angBot >= 0.0)
278 : {
279 0 : if (t2 < 0.0)
280 0 : value = 0.0;
281 0 : else if (t2 < eps)
282 0 : value = 0.5 - 0.5 * std::cos(pi * t2 / eps);
283 : else
284 0 : value = 1.0;
285 : }
286 0 : else if (angTop > 0.0 && angBot > 0.0)
287 : {
288 0 : Set::Scalar t3 = (x - nTopRight).lpNorm<2>();
289 0 : if (t3 < 0.0)
290 0 : value = 0.0;
291 0 : else if (t3 < eps)
292 0 : value = 0.5 - 0.5 * std::cos(pi * t3 / eps);
293 : else
294 0 : value = 1.0;
295 0 : }
296 : else
297 : {
298 0 : Set::Scalar t3 = (x - nBotRight).lpNorm<2>();
299 0 : if (t3 < 0.0)
300 0 : value = 0.0;
301 0 : else if (t3 < eps)
302 0 : value = 0.5 - 0.5 * std::cos(pi * t3 / eps);
303 : else
304 0 : value = 1.0;
305 : }
306 : }
307 :
308 : else
309 : {
310 0 : if (std::abs(cosRight) > cosR)
311 : {
312 0 : if (t < 0.0)
313 0 : value = 0.0;
314 0 : else if (t < eps)
315 0 : value = 0.5 - 0.5 * std::cos(pi * t / eps);
316 : else
317 0 : value = 1.0;
318 : }
319 :
320 : else
321 : {
322 0 : if (sinRight >= 0)
323 : {
324 0 : Set::Scalar t2 = (x - nTopRight).lpNorm<2>();
325 0 : if (t2 < 0.0)
326 0 : value = 0.0;
327 0 : else if (t2 < eps)
328 0 : value = 0.5 - 0.5 * std::cos(pi * t2 / eps);
329 : else
330 0 : value = 1.0;
331 : }
332 : else
333 : {
334 0 : Set::Scalar t2 = (x - nBotRight).lpNorm<2>();
335 0 : if (t2 < 0.0)
336 0 : value = 0.0;
337 0 : else if (t2 < eps)
338 0 : value = 0.5 - 0.5 * std::cos(pi * t2 / eps);
339 : else
340 0 : value = 1.0;
341 : }
342 : }
343 : }
344 : }
345 : }
346 : }
347 0 : min_value = value < min_value ? value : min_value;
348 : }
349 0 : field(i, j, k) = min_value;
350 :
351 0 : if (field(i, j, k) < 0.)
352 0 : field(i, j, k) = 0.;
353 0 : if (field(i, j, k) > 1.)
354 0 : field(i, j, k) = 1.;
355 0 : });
356 0 : }
357 0 : a_field[lev]->FillBoundary();
358 0 : }
359 :
360 : private:
361 : amrex::Vector<Set::Vector> nCenter, nOrientation, nNormal;
362 : Set::Scalar eps = 1.e-5, eps_sq = eps * eps;
363 : amrex::Vector<Set::Scalar> nThickness, nLength, nRadius;
364 : Mollifier moll = Mollifier::Dirac;
365 :
366 : public:
367 : static void
368 0 : Parse(Notch &value, IO::ParmParse &pp)
369 : {
370 0 : value.nCenter.clear();
371 0 : value.nOrientation.clear();
372 0 : value.nNormal.clear();
373 0 : value.nRadius.clear();
374 0 : value.nThickness.clear();
375 0 : value.nLength.clear();
376 :
377 0 : amrex::Vector<Set::Scalar> center;
378 0 : pp_queryarr("center", center); // Center of notch
379 :
380 0 : if (center.size() >= AMREX_SPACEDIM)
381 : {
382 0 : for (int i = 0; i < center.size(); i += AMREX_SPACEDIM)
383 0 : value.nCenter.push_back(Set::Vector(AMREX_D_DECL(center[i], center[i + 1], center[i + 2])));
384 : // AMREX_D_TERM(value.nCenter(0) = center[0];,value.nCenter(1) = center[1];,value.nCenter(2) = center[2];);
385 : }
386 : else
387 0 : Util::Abort(INFO, "Insufficient values in center");
388 :
389 0 : amrex::Vector<Set::Scalar> orientation;
390 0 : pp_queryarr("orientation", orientation); // Vector describing notch orientation
391 :
392 0 : if (orientation.size() >= AMREX_SPACEDIM && orientation.size() == center.size())
393 : {
394 0 : for (int i = 0; i < orientation.size(); i += AMREX_SPACEDIM)
395 0 : value.nOrientation.push_back(Set::Vector(AMREX_D_DECL(orientation[i], orientation[i + 1], orientation[i + 2])));
396 : // AMREX_D_TERM(value.nOrientation(0) = orientation[0];,value.nOrientation(1) = orientation[1];,value.nOrientation(2) = orientation[2];);
397 : }
398 : else
399 0 : Util::Abort(INFO, "Insufficient values in orientation");
400 :
401 0 : for (int i = 0; i < value.nOrientation.size(); i++)
402 0 : if (value.nOrientation[i].lpNorm<2>() <= 0.)
403 0 : value.nOrientation[i] = Set::Vector::Random();
404 :
405 0 : pp_queryarr("thickness", value.nThickness); // Thickness of notch
406 0 : pp_queryarr("length", value.nLength); // Length of notch
407 0 : pp_queryarr("radius", value.nRadius); // Radius of notch ends
408 0 : pp_query("eps", value.eps); // Magnitude of mollifier
409 :
410 0 : if (value.nThickness.size() == 0)
411 : {
412 0 : value.nThickness.resize(value.nCenter.size());
413 0 : for (int i = 0; i < value.nThickness.size(); i++)
414 0 : value.nThickness[i] = 0.0;
415 : }
416 0 : else if (value.nThickness.size() == 1)
417 : {
418 0 : Set::Scalar temp_thickness = value.nThickness[0] < 0. ? 0.0 : value.nThickness[0];
419 0 : value.nThickness.clear();
420 0 : value.nThickness.resize(value.nCenter.size());
421 0 : for (int i = 0; i < value.nThickness.size(); i++)
422 0 : value.nThickness[i] = temp_thickness;
423 : }
424 0 : else if (value.nThickness.size() != value.nCenter.size())
425 0 : Util::Abort(INFO, "Inconsistent size of thickness and centers");
426 : else
427 : {
428 0 : for (int i = 0; i < value.nThickness.size(); i++)
429 0 : if (value.nThickness[i] < 0.)
430 0 : value.nThickness[i] = 0.0;
431 : }
432 :
433 0 : if (value.nLength.size() != value.nCenter.size())
434 0 : Util::Abort(INFO, "Inconsistent size of length and centers");
435 0 : for (int i = 0; i < value.nLength.size(); i++)
436 0 : if (value.nLength[i] <= 0.)
437 0 : value.nLength[i] = 0.1;
438 :
439 0 : if (value.nRadius.size() == 0)
440 : {
441 0 : value.nRadius.resize(value.nThickness.size());
442 0 : for (int i = 0; i < value.nRadius.size(); i++)
443 0 : value.nRadius[i] = value.nThickness[i] / 2.0;
444 : }
445 0 : else if (value.nRadius.size() != value.nCenter.size())
446 0 : Util::Abort(INFO, "Inconsistent size of radius and centers");
447 : else
448 : {
449 0 : for (int i = 0; i < value.nRadius.size(); i++)
450 : {
451 0 : if (value.nRadius[i] <= 0.0)
452 0 : value.nRadius[i] = 0.0;
453 0 : else if (value.nRadius[i] <= value.nThickness[i])
454 0 : value.nRadius[i] = value.nThickness[i];
455 : }
456 : }
457 : // Util::Message(INFO, "value.nRadius.size() = ", value.nRadius.size(), ". value.nRadius[0] = ", value.nRadius[0]);
458 0 : if (value.eps <= 0.)
459 0 : value.eps = 1.e-5;
460 :
461 0 : std::string mollifier;
462 0 : pp_query("mollifier", mollifier); // What kind of smoother to use {dirac,gauss,erf,cos}
463 0 : if (mollifier == "Dirac" || mollifier == "dirac")
464 0 : value.moll = Mollifier::Dirac;
465 0 : else if (mollifier == "gauss" || mollifier == "Gauss")
466 : {
467 0 : value.moll = Mollifier::Gaussian;
468 0 : for (int i = 0; i < value.nRadius.size(); i++)
469 0 : value.nRadius[i] = 0.5 * value.nThickness[i];
470 : }
471 0 : else if (mollifier == "linear" || mollifier == "Linear")
472 : {
473 0 : value.moll = Mollifier::Linear;
474 0 : for (int i = 0; i < value.nRadius.size(); i++)
475 0 : value.nRadius[i] = 0.5 * value.nThickness[i];
476 : }
477 0 : else if (mollifier == "erf" || mollifier == "Erf")
478 : {
479 0 : value.moll = Mollifier::Erf;
480 0 : for (int i = 0; i < value.nRadius.size(); i++)
481 0 : value.nRadius[i] = 0.5 * value.nThickness[i];
482 : }
483 0 : else if (mollifier == "cos" || mollifier == "Cosine")
484 0 : value.moll = Mollifier::Cosine;
485 : else
486 0 : value.moll = Mollifier::Dirac;
487 :
488 0 : for (int i = 0; i < value.nOrientation.size(); i++)
489 0 : value.nOrientation[i] = value.nOrientation[i] / value.nOrientation[i].lpNorm<2>();
490 :
491 0 : value.nNormal.resize(value.nOrientation.size());
492 :
493 0 : value.eps_sq = value.eps * value.eps;
494 :
495 0 : for (int i = 0; i < value.nOrientation.size(); i++)
496 : {
497 0 : value.nNormal[i] = Set::Vector::Zero();
498 0 : if (value.nOrientation[i](0) != 0.)
499 : {
500 0 : AMREX_D_TERM(value.nNormal[i](0) = 1.;, value.nNormal[i](1) = 1.;, value.nNormal[i](2) = 1.;);
501 0 : value.nNormal[i](0) = -(AMREX_D_TERM(0., +value.nOrientation[i](1), +value.nOrientation[i](2))) / value.nOrientation[i](0);
502 0 : value.nNormal[i] = value.nNormal[i] / value.nNormal[i].lpNorm<2>();
503 : }
504 0 : else if (value.nOrientation[i](1) != 0.)
505 : {
506 0 : AMREX_D_TERM(value.nNormal[i](0) = 1.;, value.nNormal[i](1) = 1.;, value.nNormal[i](2) = 1.;);
507 0 : value.nNormal[i](1) = -(AMREX_D_TERM(value.nOrientation[i](0), +0.0, +value.nOrientation[i](2))) / value.nOrientation[i](1);
508 0 : value.nNormal[i] = value.nNormal[i] / value.nNormal[i].lpNorm<2>();
509 : }
510 : // Util::Message(INFO,"nOrientation = (", nNormal(0),",",nNormal(1),")");
511 : }
512 0 : }
513 : };
514 : }
515 : #endif
|
