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 "Set/Set.H"
12 #include "IC/IC.H"
13 
14 // Note: right now this is meant for 2D. We need to rethink this implementation for 3D.
15 namespace IC
16 {
17 class Notch : public IC
18 {
19 public:
20  enum Mollifier {Dirac, Gaussian, Cosine};
21 
22  Notch (amrex::Vector<amrex::Geometry> &_geom) : IC(_geom)
23  {
24  nCenter.resize(1); nCenter[0] = Set::Vector::Zero();
25  nOrientation.resize(1); nOrientation[0] = Set::Vector::Random();
26  nThickness.resize(1); nThickness[0] = 0.01;
27  nLength.resize(1); nLength[0] = 0.1;
28  moll = Mollifier::Dirac;
29  eps = 1e-5;
30  }
31 
32  void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
33  {
34  Set::Vector DX(geom[lev].CellSize());
35  Set::Scalar pi = std::atan(1.0)*4;
36 
37  for (amrex::MFIter mfi(*a_field[lev],true); mfi.isValid(); ++mfi)
38  {
39  //amrex::Box bx = mfi.grownnodaltilebox();
40  amrex::Box bx = mfi.tilebox();
41  bx.grow(a_field[lev]->nGrow());
42  amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
43  amrex::IndexType type = a_field[lev]->ixType();
44  amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
45 
46  Set::Vector x;
47  // NODE
48  if (type == amrex::IndexType::TheNodeType())
49  {
50  AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i)) * geom[lev].CellSize()[0];,
51  x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j)) * geom[lev].CellSize()[1];,
52  x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k)) * geom[lev].CellSize()[2];);
53  }
54  else if (type == amrex::IndexType::TheCellType())
55  {
56  AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
57  x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
58  x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
59  }
60 
61  Set::Scalar min_value = field(i,j,k);
62 
63  for (int m = 0; m < nCenter.size(); m++)
64  {
65  Set::Scalar value = 1.0;
66  if(std::abs((x-nCenter[m]).transpose()*nOrientation[m]) <= 0.5*nLength[m])
67  {
68  Set::Scalar t = std::abs((x-nCenter[m]).transpose()*nNormal[m]) - (nThickness[m]/2.0);
69  if (moll == Mollifier::Gaussian)
70  value = 0.5 + 0.5*std::erf(t / eps);
71  // value = std::erf(t/eps);
72  else if (moll == Mollifier::Dirac)
73  {
74  if (t < 0.0) value = 0.0;
75  else value = 1.0;
76  }
77  else if (moll == Mollifier::Cosine)
78  {
79  if (t < 0.0) value = 0.0;
80  else if (t < eps) value = 0.5 - 0.5*std::cos(pi*t/eps);
81  else value = 1.0;
82  }
83  }
84  else
85  {
86  Set::Vector nLeft = nCenter[m] - 0.5*nLength[m]*nOrientation[m];
87  Set::Vector nRight = nCenter[m] + 0.5*nLength[m]*nOrientation[m];
88  Set::Vector nLeft2 = nLeft, nRight2 = nRight;
89  Set::Vector correction = nOrientation[m]*(std::sqrt(nRadius[m]*nRadius[m] - nThickness[m]*nThickness[m]/4.0));
90  if (nRadius[m] > 0.5*nThickness[m])
91  {
92  nLeft2 = nLeft + correction;
93  nRight2 = nRight - correction;
94  }
95 
96  Set::Scalar cosR = (nRadius[m] <= 0.0) ? 1.0 : std::sqrt(nRadius[m]*nRadius[m] - nThickness[m]*nThickness[m]/4.0)/nRadius[m];
97  Set::Scalar cosLeft = (x-nLeft2).dot(nOrientation[m]) / ((x-nLeft2).lpNorm<2>());
98  Set::Scalar sinLeft = (x-nLeft2).dot(nNormal[m]) / ((x-nLeft2).lpNorm<2>());
99  Set::Scalar cosRight = (x-nRight2).dot(nOrientation[m]) / ((x-nRight2).lpNorm<2>());
100  Set::Scalar sinRight = (x-nRight2).dot(nNormal[m]) / ((x-nRight2).lpNorm<2>());
101 
102  Set::Scalar distLeft = (x-nLeft).dot(nOrientation[m]);
103  Set::Scalar distRight = (x-nRight).dot(nOrientation[m]);
104 
105 
106  if(distLeft <= 0.)
107  {
108  Set::Scalar t = (x-nLeft2).lpNorm<2>() - nRadius[m];
109  if (moll == Mollifier::Gaussian)
110  {
111  value = 0.5 + 0.5*std::erf( ((x-nLeft).lpNorm<2>() - nThickness[m]/2.0) / eps );
112  }
113  else if (moll == Mollifier::Dirac)
114  {
115  if (t < 0.0) value = 0.0;
116  else value = 1.0;
117  }
118  else if (moll == Mollifier::Cosine)
119  {
120  Set::Vector nTopLeft = nLeft + 0.5*nThickness[m]*nNormal[m];
121  Set::Vector nBotLeft = nLeft - 0.5*nThickness[m]*nNormal[m];
122 
123  if (nRadius[m] <= 0.0)
124  {
125  Set::Scalar t2 = -(x-nTopLeft).dot(nOrientation[m]);
126 
127  Set::Scalar angTop = (x-nTopLeft).dot(nNormal[m]);
128  Set::Scalar angBot = (x-nBotLeft).dot(nNormal[m]);
129 
130  if(angTop <= 0.0 && angBot >= 0.0)
131  {
132  if (t2 < 0.0) value = 0.0;
133  else if (t2 < eps) value = 0.5 - 0.5*std::cos(pi*t2/eps);
134  else value = 1.0;
135  }
136  else if(angTop > 0.0 && angBot > 0.0)
137  {
138  Set::Scalar t3 = (x-nTopLeft).lpNorm<2>();
139  if (t3 < 0.0) value = 0.0;
140  else if (t3 < eps) value = 0.5 - 0.5*std::cos(pi*t3/eps);
141  else value = 1.0;
142  }
143  else
144  {
145  Set::Scalar t3 = (x-nBotLeft).lpNorm<2>();
146  if (t3 < 0.0) value = 0.0;
147  else if (t3 < eps) value = 0.5 - 0.5*std::cos(pi*t3/eps);
148  else value = 1.0;
149  }
150  }
151 
152  else
153  {
154  if (std::abs(cosLeft) > cosR)
155  {
156  if (t < 0.0) value = 0.0;
157  else if (t < eps) value = 0.5 - 0.5*std::cos(pi*t/eps);
158  else value = 1.0;
159  }
160 
161  else
162  {
163  if (sinLeft >= 0)
164  {
165  Set::Scalar t2 = (x-nTopLeft).lpNorm<2>();
166  if (t2 < 0.0) value = 0.0;
167  else if (t2 < eps) value = 0.5 - 0.5*std::cos(pi*t2/eps);
168  else value = 1.0;
169  }
170  else
171  {
172  Set::Scalar t2 = (x-nBotLeft).lpNorm<2>();
173  if (t2 < 0.0) value = 0.0;
174  else if (t2 < eps) value = 0.5 - 0.5*std::cos(pi*t2/eps);
175  else value = 1.0;
176  }
177  }
178  }
179  }
180  }
181  if(distRight >= 0.)
182  {
183  Set::Scalar t = (x-nRight2).lpNorm<2>() - nRadius[m];
184  if (moll == Mollifier::Gaussian)
185  {
186  value = 0.5 + 0.5*std::erf( ((x-nRight).lpNorm<2>() - nThickness[m]/2.0) / eps );
187  }
188  else if (moll == Mollifier::Dirac)
189  {
190  if (t < 0.0) value = 0.0;
191  else value = 1.0;
192  }
193  else if (moll == Mollifier::Cosine)
194  {
195  Set::Vector nTopRight = nRight + 0.5*nThickness[m]*nNormal[m];
196  Set::Vector nBotRight = nRight - 0.5*nThickness[m]*nNormal[m];
197 
198  if (nRadius[m] <= 0.0)
199  {
200  Set::Scalar t2 = (x-nTopRight).dot(nOrientation[m]);
201  Set::Scalar angTop = (x-nTopRight).dot(nNormal[m]);
202  Set::Scalar angBot = (x-nBotRight).dot(nNormal[m]);
203 
204  if(angTop <= 0.0 && angBot >= 0.0)
205  {
206  if (t2 < 0.0) value = 0.0;
207  else if (t2 < eps) value = 0.5 - 0.5*std::cos(pi*t2/eps);
208  else value = 1.0;
209  }
210  else if(angTop > 0.0 && angBot > 0.0)
211  {
212  Set::Scalar t3 = (x-nTopRight).lpNorm<2>();
213  if (t3 < 0.0) value = 0.0;
214  else if (t3 < eps) value = 0.5 - 0.5*std::cos(pi*t3/eps);
215  else value = 1.0;
216  }
217  else
218  {
219  Set::Scalar t3 = (x-nBotRight).lpNorm<2>();
220  if (t3 < 0.0) value = 0.0;
221  else if (t3 < eps) value = 0.5 - 0.5*std::cos(pi*t3/eps);
222  else value = 1.0;
223  }
224  }
225 
226  else
227  {
228  if (std::abs(cosRight) > cosR)
229  {
230  if (t < 0.0) value = 0.0;
231  else if (t < eps) value = 0.5 - 0.5*std::cos(pi*t/eps);
232  else value = 1.0;
233  }
234 
235  else
236  {
237  if (sinRight >= 0)
238  {
239  Set::Scalar t2 = (x-nTopRight).lpNorm<2>();
240  if (t2 < 0.0) value = 0.0;
241  else if (t2 < eps) value = 0.5 - 0.5*std::cos(pi*t2/eps);
242  else value = 1.0;
243  }
244  else
245  {
246  Set::Scalar t2 = (x-nBotRight).lpNorm<2>();
247  if (t2 < 0.0) value = 0.0;
248  else if (t2 < eps) value = 0.5 - 0.5*std::cos(pi*t2/eps);
249  else value = 1.0;
250  }
251  }
252  }
253  }
254  }
255  }
256  min_value = value < min_value ? value : min_value;
257  }
258  field(i,j,k) = min_value;
259 
260  if (field(i,j,k) < 0.) field(i,j,k) = 0.;
261  if (field(i,j,k) > 1.) field(i,j,k) = 1.;
262  });
263  }
264  a_field[lev]->FillBoundary();
265  }
266 
267 private:
268  amrex::Vector<Set::Vector> nCenter, nOrientation, nNormal;
269  Set::Scalar eps = 1.e-5, eps_sq = eps*eps;
270  amrex::Vector<Set::Scalar> nThickness, nLength, nRadius;
271  Mollifier moll = Mollifier::Dirac;
272 
273 public:
274  static void Parse(Notch & value, IO::ParmParse & pp)
275  {
276  value.nCenter.clear();
277  value.nOrientation.clear();
278  value.nNormal.clear();
279  value.nRadius.clear();
280  value.nThickness.clear();
281  value.nLength.clear();
282 
283  amrex::Vector<Set::Scalar> center;
284  pp_queryarr("center",center); // Center of notch
285 
286  if(center.size() >= AMREX_SPACEDIM)
287  {
288  for (int i = 0; i<center.size(); i+=AMREX_SPACEDIM)
289  value.nCenter.push_back(Set::Vector(AMREX_D_DECL(center[i], center[i+1], center[i+2])));
290  //AMREX_D_TERM(value.nCenter(0) = center[0];,value.nCenter(1) = center[1];,value.nCenter(2) = center[2];);
291  }
292  else
293  Util::Abort(INFO, "Insufficient values in center");
294 
295  amrex::Vector<Set::Scalar> orientation;
296  pp_queryarr("orientation",orientation); // Vector describing notch orientation
297 
298  if(orientation.size()>=AMREX_SPACEDIM && orientation.size() == center.size())
299  {
300  for (int i=0; i < orientation.size(); i+=AMREX_SPACEDIM)
301  value.nOrientation.push_back(Set::Vector(AMREX_D_DECL(orientation[i], orientation[i+1], orientation[i+2])));
302  //AMREX_D_TERM(value.nOrientation(0) = orientation[0];,value.nOrientation(1) = orientation[1];,value.nOrientation(2) = orientation[2];);
303  }
304  else
305  Util::Abort(INFO, "Insufficient values in orientation");
306 
307  for (int i =0; i < value.nOrientation.size(); i++)
308  if(value.nOrientation[i].lpNorm<2>() <= 0.) value.nOrientation[i] = Set::Vector::Random();
309 
310  pp_queryarr("thickness", value.nThickness); // Thickness of notch
311  pp_queryarr("length", value.nLength); // Length of notch
312  pp_queryarr("radius", value.nRadius); // Radius of notch ends
313  pp_query("eps", value.eps); // Magnitude of mollifier
314 
315  if(value.nThickness.size() == 0)
316  {
317  value.nThickness.resize(value.nCenter.size());
318  for (int i=0; i<value.nThickness.size(); i++) value.nThickness[i] = 0.0;
319  }
320  else if(value.nThickness.size() != value.nCenter.size()) Util::Abort(INFO, "Inconsistent size of thickness and centers");
321  else
322  {
323  for (int i=0; i<value.nThickness.size(); i++)
324  if(value.nThickness[i] <= 0.) value.nThickness[i] = 0.01;
325  }
326 
327  if(value.nLength.size() != value.nCenter.size()) Util::Abort(INFO, "Inconsistent size of length and centers");
328  for (int i=0; i<value.nLength.size(); i++)
329  if(value.nLength[i] <= 0.) value.nLength[i] = 0.1;
330 
331  if (value.nRadius.size() == 0)
332  {
333  value.nRadius.resize(value.nThickness.size());
334  for (int i=0; i<value.nRadius.size(); i++)
335  value.nRadius[i] = value.nThickness[i]/2.0;
336  }
337  else if(value.nRadius.size() != value.nCenter.size()) Util::Abort(INFO, "Inconsistent size of radius and centers");
338  else
339  {
340  for (int i=0; i<value.nRadius.size(); i++)
341  {
342  if(value.nRadius[i] <= 0.0) value.nRadius[i] = 0.0;
343  else if(value.nRadius[i] <= value.nThickness[i]) value.nRadius[i] = value.nThickness[i];
344  }
345  }
346  //Util::Message(INFO, "value.nRadius.size() = ", value.nRadius.size(), ". value.nRadius[0] = ", value.nRadius[0]);
347  if(value.eps <= 0.) value.eps = 1.e-5;
348 
349  std::string mollifier;
350  pp_query("mollifier",mollifier); // What kind of smoother to use {dirac,gauss,erf,cos}
351  if(mollifier == "Dirac" || mollifier == "dirac")
352  value.moll = Mollifier::Dirac;
353  else if (mollifier == "gauss" || mollifier == "Gauss")
354  {
355  value.moll = Mollifier::Gaussian;
356  for (int i=0; i<value.nRadius.size(); i++)
357  value.nRadius[i] = 0.5*value.nThickness[i];
358  }
359  else if (mollifier == "erf" || mollifier == "Erf")
360  {
361  value.moll = Mollifier::Gaussian;
362  for (int i=0; i<value.nRadius.size(); i++)
363  value.nRadius[i] = 0.5*value.nThickness[i];
364  }
365  else if (mollifier == "cos" || mollifier == "Cosine")
366  value.moll = Mollifier::Cosine;
367  else
368  value.moll = Mollifier::Dirac;
369 
370  for (int i=0; i<value.nOrientation.size(); i++)
371  value.nOrientation[i] = value.nOrientation[i] / value.nOrientation[i].lpNorm<2>();
372 
373  value.nNormal.resize(value.nOrientation.size());
374 
375  value.eps_sq = value.eps*value.eps;
376 
377  for (int i = 0; i<value.nOrientation.size(); i++)
378  {
379  value.nNormal[i] = Set::Vector::Zero();
380  if(value.nOrientation[i](0) != 0.)
381  {
382  AMREX_D_TERM(value.nNormal[i](0) = 1.;, value.nNormal[i](1) = 1.;, value.nNormal[i](2) = 1.;);
383  value.nNormal[i](0) = -(AMREX_D_TERM(0.,+value.nOrientation[i](1),+value.nOrientation[i](2)))/value.nOrientation[i](0);
384  value.nNormal[i] = value.nNormal[i]/value.nNormal[i].lpNorm<2>();
385  }
386  else if(value.nOrientation[i](1) != 0.)
387  {
388  AMREX_D_TERM(value.nNormal[i](0) = 1.;, value.nNormal[i](1) = 1.;, value.nNormal[i](2) = 1.;);
389  value.nNormal[i](1) = -(AMREX_D_TERM(value.nOrientation[i](0), + 0.0, + value.nOrientation[i](2)))/value.nOrientation[i](1);
390  value.nNormal[i] = value.nNormal[i]/value.nNormal[i].lpNorm<2>();
391  }
392  //Util::Message(INFO,"nOrientation = (", nNormal(0),",",nNormal(1),")");
393  }
394  }
395 };
396 }
397 #endif
IC::IC::geom
amrex::Vector< amrex::Geometry > & geom
Definition: IC.H:45
IC::Notch::Notch
Notch(amrex::Vector< amrex::Geometry > &_geom)
Definition: Notch.H:22
BC::AMREX_D_DECL
@ AMREX_D_DECL
Definition: BC.H:34
t
std::time_t t
Definition: FileNameParse.cpp:12
Set::Field< Set::Scalar >
Definition: Set.H:236
IC::Notch::nNormal
amrex::Vector< Set::Vector > nNormal
Definition: Notch.H:268
IC::Notch::Cosine
@ Cosine
Definition: Notch.H:20
IC::Notch::Dirac
@ Dirac
Definition: Notch.H:20
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition: Base.H:20
pp_query
#define pp_query(...)
Definition: ParmParse.H:104
Util::Gaussian
Set::Scalar Gaussian(amrex::Real mean, amrex::Real std_deviation)
Definition: Set.cpp:13
Util::Random
Set::Scalar Random()
Definition: Set.cpp:9
Set::Scalar
amrex::Real Scalar
Definition: Base.H:19
IC::Notch
Definition: Notch.H:17
IC::Notch::nThickness
amrex::Vector< Set::Scalar > nThickness
Definition: Notch.H:270
IC::Notch::nCenter
amrex::Vector< Set::Vector > nCenter
Definition: Notch.H:268
IC::Notch::nLength
amrex::Vector< Set::Scalar > nLength
Definition: Notch.H:270
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:165
IC::Notch::Parse
static void Parse(Notch &value, IO::ParmParse &pp)
Definition: Notch.H:274
IC::Notch::nRadius
amrex::Vector< Set::Scalar > nRadius
Definition: Notch.H:270
IC::Notch::eps_sq
Set::Scalar eps_sq
Definition: Notch.H:269
Set.H
IC
Definition: Affine.H:18
IC::Notch::Mollifier
Mollifier
Definition: Notch.H:20
IO::ParmParse
Definition: ParmParse.H:110
IC::Notch::Gaussian
@ Gaussian
Definition: Notch.H:20
IC::Notch::nOrientation
amrex::Vector< Set::Vector > nOrientation
Definition: Notch.H:268
IC::Notch::eps
Set::Scalar eps
Definition: Notch.H:269
IC::Notch::Add
void Add(const int &lev, Set::Field< Set::Scalar > &a_field, Set::Scalar)
Definition: Notch.H:32
INFO
#define INFO
Definition: Util.H:20
IC.H
IC::Notch::moll
Mollifier moll
Definition: Notch.H:271
pp_queryarr
#define pp_queryarr(...)
Definition: ParmParse.H:103
Generated on Fri Aug 16 2024 16:30:29 for Alamo by   doxygen 1.8.17