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