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 if (center.size() >= AMREX_SPACEDIM)
381 {
382 for (int i = 0; i < center.size(); i += AMREX_SPACEDIM)
383 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 Util::Abort(INFO, "Insufficient values in center");
388
389 amrex::Vector<Set::Scalar> orientation;
390 pp_queryarr("orientation", orientation); // Vector describing notch orientation
391
392 if (orientation.size() >= AMREX_SPACEDIM && orientation.size() == center.size())
393 {
394 for (int i = 0; i < orientation.size(); i += AMREX_SPACEDIM)
395 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 Util::Abort(INFO, "Insufficient values in orientation");
400
401 for (int i = 0; i < value.nOrientation.size(); i++)
402 if (value.nOrientation[i].lpNorm<2>() <= 0.)
403 value.nOrientation[i] = Set::Vector::Random();
404
405 pp_queryarr("thickness", value.nThickness); // Thickness of notch
406 pp_queryarr("length", value.nLength); // Length of notch
407 pp_queryarr("radius", value.nRadius); // Radius of notch ends
408 pp_query("eps", value.eps); // Magnitude of mollifier
409
410 if (value.nThickness.size() == 0)
411 {
412 value.nThickness.resize(value.nCenter.size());
413 for (int i = 0; i < value.nThickness.size(); i++)
414 value.nThickness[i] = 0.0;
415 }
416 else if (value.nThickness.size() == 1)
417 {
418 Set::Scalar temp_thickness = value.nThickness[0] < 0. ? 0.0 : value.nThickness[0];
419 value.nThickness.clear();
420 value.nThickness.resize(value.nCenter.size());
421 for (int i = 0; i < value.nThickness.size(); i++)
422 value.nThickness[i] = temp_thickness;
423 }
424 else if (value.nThickness.size() != value.nCenter.size())
425 Util::Abort(INFO, "Inconsistent size of thickness and centers");
426 else
427 {
428 for (int i = 0; i < value.nThickness.size(); i++)
429 if (value.nThickness[i] < 0.)
430 value.nThickness[i] = 0.0;
431 }
432
433 if (value.nLength.size() != value.nCenter.size())
434 Util::Abort(INFO, "Inconsistent size of length and centers");
435 for (int i = 0; i < value.nLength.size(); i++)
436 if (value.nLength[i] <= 0.)
437 value.nLength[i] = 0.1;
438
439 if (value.nRadius.size() == 0)
440 {
441 value.nRadius.resize(value.nThickness.size());
442 for (int i = 0; i < value.nRadius.size(); i++)
443 value.nRadius[i] = value.nThickness[i] / 2.0;
444 }
445 else if (value.nRadius.size() != value.nCenter.size())
446 Util::Abort(INFO, "Inconsistent size of radius and centers");
447 else
448 {
449 for (int i = 0; i < value.nRadius.size(); i++)
450 {
451 if (value.nRadius[i] <= 0.0)
452 value.nRadius[i] = 0.0;
453 else if (value.nRadius[i] <= value.nThickness[i])
454 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 if (value.eps <= 0.)
459 value.eps = 1.e-5;
460
461 std::string mollifier;
462 pp_query("mollifier", mollifier); // What kind of smoother to use {dirac,gauss,erf,cos}
463 if (mollifier == "Dirac" || mollifier == "dirac")
464 value.moll = Mollifier::Dirac;
465 else if (mollifier == "gauss" || mollifier == "Gauss")
466 {
467 value.moll = Mollifier::Gaussian;
468 for (int i = 0; i < value.nRadius.size(); i++)
469 value.nRadius[i] = 0.5 * value.nThickness[i];
470 }
471 else if (mollifier == "linear" || mollifier == "Linear")
472 {
473 value.moll = Mollifier::Linear;
474 for (int i = 0; i < value.nRadius.size(); i++)
475 value.nRadius[i] = 0.5 * value.nThickness[i];
476 }
477 else if (mollifier == "erf" || mollifier == "Erf")
478 {
479 value.moll = Mollifier::Erf;
480 for (int i = 0; i < value.nRadius.size(); i++)
481 value.nRadius[i] = 0.5 * value.nThickness[i];
482 }
483 else if (mollifier == "cos" || mollifier == "Cosine")
484 value.moll = Mollifier::Cosine;
485 else
486 value.moll = Mollifier::Dirac;
487
488 for (int i = 0; i < value.nOrientation.size(); i++)
489 value.nOrientation[i] = value.nOrientation[i] / value.nOrientation[i].lpNorm<2>();
490
491 value.nNormal.resize(value.nOrientation.size());
492
493 value.eps_sq = value.eps * value.eps;
494
495 for (int i = 0; i < value.nOrientation.size(); i++)
496 {
497 value.nNormal[i] = Set::Vector::Zero();
498 if (value.nOrientation[i](0) != 0.)
499 {
500 AMREX_D_TERM(value.nNormal[i](0) = 1.;, value.nNormal[i](1) = 1.;, value.nNormal[i](2) = 1.;);
501 value.nNormal[i](0) = -(AMREX_D_TERM(0., +value.nOrientation[i](1), +value.nOrientation[i](2))) / value.nOrientation[i](0);
502 value.nNormal[i] = value.nNormal[i] / value.nNormal[i].lpNorm<2>();
503 }
504 else if (value.nOrientation[i](1) != 0.)
505 {
506 AMREX_D_TERM(value.nNormal[i](0) = 1.;, value.nNormal[i](1) = 1.;, value.nNormal[i](2) = 1.;);
507 value.nNormal[i](1) = -(AMREX_D_TERM(value.nOrientation[i](0), +0.0, +value.nOrientation[i](2))) / value.nOrientation[i](1);
508 value.nNormal[i] = value.nNormal[i] / value.nNormal[i].lpNorm<2>();
509 }
510 // Util::Message(INFO,"nOrientation = (", nNormal(0),",",nNormal(1),")");
511 }
512 }
513};
514}
515#endif
t
std::time_t t
Definition FileNameParse.cpp:9
pp_queryarr
#define pp_queryarr(...)
Definition ParmParse.H:105
pp_query
#define pp_query(...)
Definition ParmParse.H:108
pp_queryclass
#define pp_queryclass(...)
Definition ParmParse.H:109
INFO
#define INFO
Definition Util.H:21
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
IC
Initialize a spherical inclusion.
Definition BMP.H:20
Set
A collection of data types and symmetry-reduced data structures.
Definition Base.H:17
Set::Scalar
amrex::Real Scalar
Definition Base.H:18
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:19
Util::Abort
void Abort(const char *msg)
Definition Util.cpp:225
Generated by doxygen 1.9.8