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