Alamo
Expression.H
Go to the documentation of this file.
1 #ifndef BC_OPERATOR_ELASTIC_EXPRESSION_H
2 #define BC_OPERATOR_ELASTIC_EXPRESSION_H
3 
4 // #include "Operator/Elastic.H"
5 #include "IO/ParmParse.H"
6 #include "Numeric/Interpolator/Linear.H"
7 #include "BC/Operator/Elastic/Elastic.H"
8 
9 #include "AMReX_Parser.H"
10 
11 namespace BC
12 {
13 namespace Operator
14 {
15 namespace Elastic
16 {
17 class Expression : public Elastic
18 {
19 public:
20  static constexpr const char* name = "expression";
21 
22  //static constexpr const char* const Elastic::strings[];
23  #if AMREX_SPACEDIM==2
24  static const constexpr char * const facestr[] = {
25  "xlo","ylo","xhi","yhi","xloylo","xloyhi","xhiylo","xhiyhi","int"
26  };
27  #elif AMREX_SPACEDIM==3
28  static const constexpr char * const facestr[] = {
29  "xlo","ylo","zlo","xhi","yhi","zhi",
30  "ylozlo","ylozhi","yhizlo","yhizhi",
31  "zloxlo","zloxhi","zhixlo","zhixhi",
32  "xloylo","xloyhi","xhiylo","xhiyhi",
33  "xloylozlo","xloylozhi","xloyhizlo","xloyhizhi",
34  "xhiylozlo","xhiylozhi","xhiyhizlo","xhiyhizhi",
35  "int"
36  };
37  #endif
38 
39  Expression()
40  {
41  // By default, all boundary conditions are displacement
42  for (int face = 0; face < m_nfaces; face++)
43  for (int direction = 0; direction < AMREX_SPACEDIM; direction++)
44  {
45  m_bc_type[face][direction] = Type::Displacement;
46  }
47  };
48  ~Expression() {};
49 
50 
51  using Elastic::Init;
52 
53  virtual void
54  Init(amrex::FabArray<amrex::BaseFab<Set::Vector>> *a_rhs,
55  const amrex::Geometry &a_geom,
56  bool a_homogeneous = false) const override
57  {
58  amrex::Box domain(a_geom.Domain());
59  domain.convert(amrex::IntVect::TheNodeVector());
60  const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
61  for (amrex::MFIter mfi(*a_rhs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
62  {
63  amrex::Box bx = mfi.tilebox();
64  bx.grow(2);
65  bx = bx & domain;
66  amrex::Array4<Set::Vector> const& rhs = a_rhs->array(mfi);
67  amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
68 
69  for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
70  {
71  Face face = Face::INT;
72 
73  #if AMREX_SPACEDIM == 2
74 
75  if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
76  else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
77  else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
78  else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
79 
80  else if (i==lo.x) face = Face::XLO;
81  else if (i==hi.x) face = Face::XHI;
82  else if (j==lo.y) face = Face::YLO;
83  else if (j==hi.y) face = Face::YHI;
84 
85  #elif AMREX_SPACEDIM == 3
86 
87  if (i==lo.x && j==lo.y && k==lo.z) face = Face::XLO_YLO_ZLO;
88  else if (i==lo.x && j==lo.y && k==hi.z) face = Face::XLO_YLO_ZHI;
89  else if (i==lo.x && j==hi.y && k==lo.z) face = Face::XLO_YHI_ZLO;
90  else if (i==lo.x && j==hi.y && k==hi.z) face = Face::XLO_YHI_ZHI;
91  else if (i==hi.x && j==lo.y && k==lo.z) face = Face::XHI_YLO_ZLO;
92  else if (i==hi.x && j==lo.y && k==hi.z) face = Face::XHI_YLO_ZHI;
93  else if (i==hi.x && j==hi.y && k==lo.z) face = Face::XHI_YHI_ZLO;
94  else if (i==hi.x && j==hi.y && k==hi.z) face = Face::XHI_YHI_ZHI;
95 
96  else if (j==lo.y && k==lo.z) face = Face::YLO_ZLO;
97  else if (j==lo.y && k==hi.z) face = Face::YLO_ZHI;
98  else if (j==hi.y && k==lo.z) face = Face::YHI_ZLO;
99  else if (j==hi.y && k==hi.z) face = Face::YHI_ZHI;
100  else if (k==lo.z && i==lo.x) face = Face::ZLO_XLO;
101  else if (k==lo.z && i==hi.x) face = Face::ZLO_XHI;
102  else if (k==hi.z && i==lo.x) face = Face::ZHI_XLO;
103  else if (k==hi.z && i==hi.x) face = Face::ZHI_XHI;
104  else if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
105  else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
106  else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
107  else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
108 
109  else if (i==lo.x) face = Face::XLO;
110  else if (i==hi.x) face = Face::XHI;
111  else if (j==lo.y) face = Face::YLO;
112  else if (j==hi.y) face = Face::YHI;
113  else if (k==lo.z) face = Face::ZLO;
114  else if (k==hi.z) face = Face::ZHI;
115 
116  #endif
117 
118  amrex::IndexType type = a_rhs->ixType();
119  if (!(face == Face::INT))
120  {
121  if (a_homogeneous && m_bc_type[face][dir] == Type::Displacement)
122  rhs(i,j,k)(dir) = 0.0;
123  else
124  {
125  Set::Vector x = Set::Position(i, j, k, a_geom, type);
126  #if AMREX_SPACEDIM==1
127  rhs(i,j,k)(dir) = m_bc_func[face][dir](x(0),0.0,0.0,m_time);
128  #elif AMREX_SPACEDIM==2
129  rhs(i,j,k)(dir) = m_bc_func[face][dir](x(0),x(1),0.0,m_time);
130  #elif AMREX_SPACEDIM==3
131  rhs(i,j,k)(dir) = m_bc_func[face][dir](x(0),x(1),x(2),m_time);
132  #endif
133  }
134  }
135  }
136  });
137  }
138  }
139 
140 
141  virtual void
142  Init(amrex::MultiFab * a_rhs,
143  const amrex::Geometry &a_geom,
144  bool a_homogeneous = false) const override
145  {
146  amrex::Box domain(a_geom.Domain());
147  domain.convert(amrex::IntVect::TheNodeVector());
148  const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
149  for (amrex::MFIter mfi(*a_rhs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
150  {
151  amrex::Box bx = mfi.tilebox();
152  bx.grow(2);
153  bx = bx & domain;
154  amrex::Array4<amrex::Real> const& rhs = a_rhs->array(mfi);
155  amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
156 
157  for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
158  {
159  Face face = Face::INT;
160 
161  #if AMREX_SPACEDIM == 2
162 
163  if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
164  else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
165  else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
166  else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
167 
168  else if (i==lo.x) face = Face::XLO;
169  else if (i==hi.x) face = Face::XHI;
170  else if (j==lo.y) face = Face::YLO;
171  else if (j==hi.y) face = Face::YHI;
172 
173  #elif AMREX_SPACEDIM == 3
174 
175  if (i==lo.x && j==lo.y && k==lo.z) face = Face::XLO_YLO_ZLO;
176  else if (i==lo.x && j==lo.y && k==hi.z) face = Face::XLO_YLO_ZHI;
177  else if (i==lo.x && j==hi.y && k==lo.z) face = Face::XLO_YHI_ZLO;
178  else if (i==lo.x && j==hi.y && k==hi.z) face = Face::XLO_YHI_ZHI;
179  else if (i==hi.x && j==lo.y && k==lo.z) face = Face::XHI_YLO_ZLO;
180  else if (i==hi.x && j==lo.y && k==hi.z) face = Face::XHI_YLO_ZHI;
181  else if (i==hi.x && j==hi.y && k==lo.z) face = Face::XHI_YHI_ZLO;
182  else if (i==hi.x && j==hi.y && k==hi.z) face = Face::XHI_YHI_ZHI;
183 
184  else if (j==lo.y && k==lo.z) face = Face::YLO_ZLO;
185  else if (j==lo.y && k==hi.z) face = Face::YLO_ZHI;
186  else if (j==hi.y && k==lo.z) face = Face::YHI_ZLO;
187  else if (j==hi.y && k==hi.z) face = Face::YHI_ZHI;
188  else if (k==lo.z && i==lo.x) face = Face::ZLO_XLO;
189  else if (k==lo.z && i==hi.x) face = Face::ZLO_XHI;
190  else if (k==hi.z && i==lo.x) face = Face::ZHI_XLO;
191  else if (k==hi.z && i==hi.x) face = Face::ZHI_XHI;
192  else if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
193  else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
194  else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
195  else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
196 
197  else if (i==lo.x) face = Face::XLO;
198  else if (i==hi.x) face = Face::XHI;
199  else if (j==lo.y) face = Face::YLO;
200  else if (j==hi.y) face = Face::YHI;
201  else if (k==lo.z) face = Face::ZLO;
202  else if (k==hi.z) face = Face::ZHI;
203 
204  #endif
205 
206  amrex::IndexType type = a_rhs->ixType();
207  if (!(face == Face::INT))
208  {
209  if (a_homogeneous && m_bc_type[face][dir] == Type::Displacement)
210  rhs(i,j,k,dir) = 0.0;
211  else
212  {
213  Set::Vector x = Set::Position(i, j, k, a_geom, type);
214  #if AMREX_SPACEDIM==1
215  rhs(i,j,k,dir) = (m_bc_func[face][dir])(x(0),0.0,0.0,m_time);
216  #elif AMREX_SPACEDIM==2
217  rhs(i,j,k,dir) = (m_bc_func[face][dir])(x(0),x(1),0.0,m_time);
218  #elif AMREX_SPACEDIM==3
219  rhs(i,j,k,dir) = (m_bc_func[face][dir])(x(0),x(1),x(2),m_time);
220  #endif
221  }
222  }
223  }
224  });
225  }
226  }
227 
228 
229 
230  virtual
231  std::array<Type,AMREX_SPACEDIM> getType (
232  const int &i, const int &j, [[maybe_unused]] const int &k,
233  const amrex::Box &domain) override
234  {
235  (void)i, (void)j, (void)k;
236  amrex::IntVect m(AMREX_D_DECL(i,j,k));
237  const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
238 
239  // Corners
240  #if AMREX_SPACEDIM == 2
241  if (m[0] == lo.x && m[1] == lo.y) return m_bc_type[Face::XLO_YLO];
242  if (m[0] == lo.x && m[1] == hi.y) return m_bc_type[Face::XLO_YHI];
243  if (m[0] == hi.x && m[1] == lo.y) return m_bc_type[Face::XHI_YLO];
244  if (m[0] == hi.x && m[1] == hi.y) return m_bc_type[Face::XHI_YHI];
245  if (m[0] == lo.x) return m_bc_type[Face::XLO];
246  if (m[0] == hi.x) return m_bc_type[Face::XHI];
247  if (m[1] == lo.y) return m_bc_type[Face::YLO];
248  if (m[1] == hi.y) return m_bc_type[Face::YHI];
249 
250  #elif AMREX_SPACEDIM == 3
251 
252  if (m[0] == lo.x && m[1] == lo.y && m[2] == lo.z) return m_bc_type[Face::XLO_YLO_ZLO];
253  if (m[0] == lo.x && m[1] == lo.y && m[2] == hi.z) return m_bc_type[Face::XLO_YLO_ZHI];
254  if (m[0] == lo.x && m[1] == hi.y && m[2] == lo.z) return m_bc_type[Face::XLO_YHI_ZLO];
255  if (m[0] == lo.x && m[1] == hi.y && m[2] == hi.z) return m_bc_type[Face::XLO_YHI_ZHI];
256  if (m[0] == hi.x && m[1] == lo.y && m[2] == lo.z) return m_bc_type[Face::XHI_YLO_ZLO];
257  if (m[0] == hi.x && m[1] == lo.y && m[2] == hi.z) return m_bc_type[Face::XHI_YLO_ZHI];
258  if (m[0] == hi.x && m[1] == hi.y && m[2] == lo.z) return m_bc_type[Face::XHI_YHI_ZLO];
259  if (m[0] == hi.x && m[1] == hi.y && m[2] == hi.z) return m_bc_type[Face::XHI_YHI_ZHI];
260 
261  if (m[1] == lo.y && m[2] == lo.z) return m_bc_type[Face::YLO_ZLO];
262  if (m[1] == lo.y && m[2] == hi.z) return m_bc_type[Face::YLO_ZHI];
263  if (m[1] == hi.y && m[2] == lo.z) return m_bc_type[Face::YHI_ZLO];
264  if (m[1] == hi.y && m[2] == hi.z) return m_bc_type[Face::YHI_ZHI];
265  if (m[2] == lo.z && m[0] == lo.x) return m_bc_type[Face::ZLO_XLO];
266  if (m[2] == lo.z && m[0] == hi.x) return m_bc_type[Face::ZLO_XHI];
267  if (m[2] == hi.z && m[0] == lo.x) return m_bc_type[Face::ZHI_XLO];
268  if (m[2] == hi.z && m[0] == hi.x) return m_bc_type[Face::ZHI_XHI];
269  if (m[0] == lo.x && m[1] == lo.y) return m_bc_type[Face::XLO_YLO];
270  if (m[0] == lo.x && m[1] == hi.y) return m_bc_type[Face::XLO_YHI];
271  if (m[0] == hi.x && m[1] == lo.y) return m_bc_type[Face::XHI_YLO];
272  if (m[0] == hi.x && m[1] == hi.y) return m_bc_type[Face::XHI_YHI];
273 
274  if (m[0] == lo.x) return m_bc_type[Face::XLO];
275  if (m[0] == hi.x) return m_bc_type[Face::XHI];
276  if (m[1] == lo.y) return m_bc_type[Face::YLO];
277  if (m[1] == hi.y) return m_bc_type[Face::YHI];
278  if (m[2] == lo.z) return m_bc_type[Face::ZLO];
279  if (m[2] == hi.z) return m_bc_type[Face::ZHI];
280 
281  #endif
282 
283  return {AMREX_D_DECL(Type::None,Type::None,Type::None)};
284  }
285 
286 
287  AMREX_FORCE_INLINE
288  Set::Vector operator () (const Set::Vector &u,
289  const Set::Matrix &gradu,
290  const Set::Matrix &sigma,
291  const int &i, const int &j, const int &k,
292  const amrex::Box &domain) override
293  {
294  (void)i; (void)j; (void)k; // Suppress "unused variable" warnings
295  //Set::Vector f;
296 
297  const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
298 
299  amrex::IntVect m(AMREX_D_DECL(i,j,k));
300 
301  // Corners
302  #if AMREX_SPACEDIM == 2
303 
304  if (m[0] == lo.x && m[1] == lo.y) return set(m_bc_type[Face::XLO_YLO], u, gradu, sigma, Set::Vector(-SQRT2INV, -SQRT2INV));
305  if (m[0] == lo.x && m[1] == hi.y) return set(m_bc_type[Face::XLO_YHI], u, gradu, sigma, Set::Vector(-SQRT2INV, +SQRT2INV));
306  if (m[0] == hi.x && m[1] == lo.y) return set(m_bc_type[Face::XHI_YLO], u, gradu, sigma, Set::Vector(+SQRT2INV, -SQRT2INV));
307  if (m[0] == hi.x && m[1] == hi.y) return set(m_bc_type[Face::XHI_YHI], u, gradu, sigma, Set::Vector(+SQRT2INV, +SQRT2INV));
308 
309  if (m[0] == lo.x) return set(m_bc_type[Face::XLO], u, gradu, sigma, Set::Vector(-1, 0));
310  if (m[0] == hi.x) return set(m_bc_type[Face::XHI], u, gradu, sigma, Set::Vector(+1, 0));
311  if (m[1] == lo.y) return set(m_bc_type[Face::YLO], u, gradu, sigma, Set::Vector( 0,-1));
312  if (m[1] == hi.y) return set(m_bc_type[Face::YHI], u, gradu, sigma, Set::Vector( 0,+1));
313 
314  #elif AMREX_SPACEDIM == 3
315 
316  if (m[0] == lo.x && m[1] == lo.y && m[2] == lo.z) return set(m_bc_type[Face::XLO_YLO_ZLO], u, gradu, sigma, Set::Vector(-SQRT3INV,-SQRT3INV,-SQRT3INV));
317  if (m[0] == lo.x && m[1] == lo.y && m[2] == hi.z) return set(m_bc_type[Face::XLO_YLO_ZHI], u, gradu, sigma, Set::Vector(-SQRT3INV,-SQRT3INV,+SQRT3INV));
318  if (m[0] == lo.x && m[1] == hi.y && m[2] == lo.z) return set(m_bc_type[Face::XLO_YHI_ZLO], u, gradu, sigma, Set::Vector(-SQRT3INV,+SQRT3INV,-SQRT3INV));
319  if (m[0] == lo.x && m[1] == hi.y && m[2] == hi.z) return set(m_bc_type[Face::XLO_YHI_ZHI], u, gradu, sigma, Set::Vector(-SQRT3INV,+SQRT3INV,+SQRT3INV));
320  if (m[0] == hi.x && m[1] == lo.y && m[2] == lo.z) return set(m_bc_type[Face::XHI_YLO_ZLO], u, gradu, sigma, Set::Vector(+SQRT3INV,-SQRT3INV,-SQRT3INV));
321  if (m[0] == hi.x && m[1] == lo.y && m[2] == hi.z) return set(m_bc_type[Face::XHI_YLO_ZHI], u, gradu, sigma, Set::Vector(+SQRT3INV,-SQRT3INV,+SQRT3INV));
322  if (m[0] == hi.x && m[1] == hi.y && m[2] == lo.z) return set(m_bc_type[Face::XHI_YHI_ZLO], u, gradu, sigma, Set::Vector(+SQRT3INV,+SQRT3INV,-SQRT3INV));
323  if (m[0] == hi.x && m[1] == hi.y && m[2] == hi.z) return set(m_bc_type[Face::XHI_YHI_ZHI], u, gradu, sigma, Set::Vector(+SQRT3INV,+SQRT3INV,+SQRT3INV));
324 
325  if (m[1] == lo.y && m[2] == lo.z) return set(m_bc_type[Face::YLO_ZLO], u, gradu, sigma, Set::Vector(0, -SQRT2INV,-SQRT2INV));
326  if (m[1] == lo.y && m[2] == hi.z) return set(m_bc_type[Face::YLO_ZHI], u, gradu, sigma, Set::Vector(0, -SQRT2INV,+SQRT2INV));
327  if (m[1] == hi.y && m[2] == lo.z) return set(m_bc_type[Face::YHI_ZLO], u, gradu, sigma, Set::Vector(0, +SQRT2INV,-SQRT2INV));
328  if (m[1] == hi.y && m[2] == hi.z) return set(m_bc_type[Face::YHI_ZHI], u, gradu, sigma, Set::Vector(0, +SQRT2INV,+SQRT2INV));
329  if (m[2] == lo.z && m[0] == lo.x) return set(m_bc_type[Face::ZLO_XLO], u, gradu, sigma, Set::Vector(-SQRT2INV, 0, -SQRT2INV));
330  if (m[2] == lo.z && m[0] == hi.x) return set(m_bc_type[Face::ZLO_XHI], u, gradu, sigma, Set::Vector(+SQRT2INV, 0, -SQRT2INV));
331  if (m[2] == hi.z && m[0] == lo.x) return set(m_bc_type[Face::ZHI_XLO], u, gradu, sigma, Set::Vector(-SQRT2INV, 0, +SQRT2INV));
332  if (m[2] == hi.z && m[0] == hi.x) return set(m_bc_type[Face::ZHI_XHI], u, gradu, sigma, Set::Vector(+SQRT2INV, 0, +SQRT2INV));
333  if (m[0] == lo.x && m[1] == lo.y) return set(m_bc_type[Face::XLO_YLO], u, gradu, sigma, Set::Vector(-SQRT2INV, -SQRT2INV,0 ));
334  if (m[0] == lo.x && m[1] == hi.y) return set(m_bc_type[Face::XLO_YHI], u, gradu, sigma, Set::Vector(-SQRT2INV, +SQRT2INV,0 ));
335  if (m[0] == hi.x && m[1] == lo.y) return set(m_bc_type[Face::XHI_YLO], u, gradu, sigma, Set::Vector(+SQRT2INV, -SQRT2INV,0 ));
336  if (m[0] == hi.x && m[1] == hi.y) return set(m_bc_type[Face::XHI_YHI], u, gradu, sigma, Set::Vector(+SQRT2INV, +SQRT2INV,0 ));
337 
338  if (m[0] == lo.x) return set(m_bc_type[Face::XLO], u, gradu, sigma, Set::Vector(-1, 0, 0));
339  if (m[0] == hi.x) return set(m_bc_type[Face::XHI], u, gradu, sigma, Set::Vector(+1, 0, 0));
340  if (m[1] == lo.y) return set(m_bc_type[Face::YLO], u, gradu, sigma, Set::Vector( 0,-1, 0));
341  if (m[1] == hi.y) return set(m_bc_type[Face::YHI], u, gradu, sigma, Set::Vector( 0,+1, 0));
342  if (m[2] == lo.z) return set(m_bc_type[Face::ZLO], u, gradu, sigma, Set::Vector( 0, 0,-1));
343  if (m[2] == hi.z) return set(m_bc_type[Face::ZHI], u, gradu, sigma, Set::Vector( 0, 0,+1));
344 
345  #endif
346 
347  Util::Abort(INFO,"Boundary condition error");
348  return Set::Vector::Zero();
349  }
350 
351  AMREX_FORCE_INLINE
352  Set::Vector set(std::array<Type,AMREX_SPACEDIM> &bc_type,
353  const Set::Vector &u, const Set::Matrix &gradu, const Set::Matrix &sigma, Set::Vector n) const
354  {
355  Set::Vector f = Set::Vector::Zero();
356  for (int i = 0; i < AMREX_SPACEDIM; i++)
357  {
358  if (bc_type[i] == Type::Displacement)
359  f(i) = u(i);
360  else if (bc_type[i] == Type::Traction)
361  f(i) = (sigma*n)(i);
362  else if (bc_type[i] == Type::Neumann)
363  f(i) = (gradu*n)(i);
364  else if (bc_type[i] == Periodic)
365  continue;
366  }
367  return f;
368  }
369 
370 protected:
371 
372  #if AMREX_SPACEDIM==2
373  static const int m_nfaces = 8;
374  #elif AMREX_SPACEDIM==3
375  static const int m_nfaces = 26;
376  #endif
377 
378  std::array<std::array<Type, AMREX_SPACEDIM>, m_nfaces> m_bc_type;
379  //std::array<std::array<FunctionParserAD,AMREX_SPACEDIM>, m_nfaces> m_bc_func;
380  std::array<std::array<amrex::Parser,AMREX_SPACEDIM>, m_nfaces> m_bc_func_parser;
381  std::array<std::array<amrex::ParserExecutor<4>,AMREX_SPACEDIM>, m_nfaces> m_bc_func;
382 
383 public:
384 
385  //
386  // This is basically the same as :ref:`BC::Operator::Elastic::Constant` except that
387  // you can use dynamically compiled expressions in space and time to define values.
388  //
389  // Usage is exactly the same except that the "val" inputs can depend on x, y, z, and t.
390  //
391  static void Parse(Expression & value, IO::ParmParse & pp)
392  {
393  std::map<std::string, Type> bcmap;
394  bcmap["displacement"] = Type::Displacement;
395  bcmap["disp"] = Type::Displacement;
396  bcmap["neumann"] = Type::Neumann;
397  bcmap["traction"] = Type::Traction;
398  bcmap["trac"] = Type::Traction;
399  bcmap["periodic"] = Type::Periodic;
400 
401  std::vector<std::string> str;
402 
403  // TYPES
404 
405  #if AMREX_SPACEDIM==3
406  if (pp.contains("type.xloylozlo")) { pp_queryarr("type.xloylozlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XLO_YLO_ZLO][i] = bcmap[str[i]]; } // 3D Corner
407  if (pp.contains("type.xloylozhi")) { pp_queryarr("type.xloylozhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XLO_YLO_ZHI][i] = bcmap[str[i]]; } // 3D Corner
408  if (pp.contains("type.xloyhizlo")) { pp_queryarr("type.xloyhizlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XLO_YHI_ZLO][i] = bcmap[str[i]]; } // 3D Corner
409  if (pp.contains("type.xloyhizhi")) { pp_queryarr("type.xloyhizhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XLO_YHI_ZHI][i] = bcmap[str[i]]; } // 3D Corner
410  if (pp.contains("type.xhiylozlo")) { pp_queryarr("type.xhiylozlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XHI_YLO_ZLO][i] = bcmap[str[i]]; } // 3D Corner
411  if (pp.contains("type.xhiylozhi")) { pp_queryarr("type.xhiylozhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XHI_YLO_ZHI][i] = bcmap[str[i]]; } // 3D Corner
412  if (pp.contains("type.xhiyhizlo")) { pp_queryarr("type.xhiyhizlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XHI_YHI_ZLO][i] = bcmap[str[i]]; } // 3D Corner
413  if (pp.contains("type.xhiyhizhi")) { pp_queryarr("type.xhiyhizhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XHI_YHI_ZHI][i] = bcmap[str[i]]; } // 3D Corner
414  #endif
415 
416  #if AMREX_SPACEDIM==3
417  if (pp.contains("type.ylozlo")) { pp_queryarr("type.ylozlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::YLO_ZLO][i] = bcmap[str[i]]; } // 3D Edge
418  if (pp.contains("type.ylozhi")) { pp_queryarr("type.ylozhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::YLO_ZHI][i] = bcmap[str[i]]; } // 3D Edge
419  if (pp.contains("type.yhizlo")) { pp_queryarr("type.yhizlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::YHI_ZLO][i] = bcmap[str[i]]; } // 3D Edge
420  if (pp.contains("type.yhizhi")) { pp_queryarr("type.yhizhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::YHI_ZHI][i] = bcmap[str[i]]; } // 3D Edge
421  if (pp.contains("type.zloxlo")) { pp_queryarr("type.zloxlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::ZLO_XLO][i] = bcmap[str[i]]; } // 3D Edge
422  if (pp.contains("type.zloxhi")) { pp_queryarr("type.zloxhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::ZLO_XHI][i] = bcmap[str[i]]; } // 3D Edge
423  if (pp.contains("type.zhixlo")) { pp_queryarr("type.zhixlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::ZHI_XLO][i] = bcmap[str[i]]; } // 3D Edge
424  if (pp.contains("type.zhixhi")) { pp_queryarr("type.zhixhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::ZHI_XHI][i] = bcmap[str[i]]; } // 3D Edge
425  #endif
426  if (pp.contains("type.xloylo")) { pp_queryarr("type.xloylo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XLO_YLO][i] = bcmap[str[i]]; } // 3D Edge / 2D Corner
427  if (pp.contains("type.xloyhi")) { pp_queryarr("type.xloyhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XLO_YHI][i] = bcmap[str[i]]; } // 3D Edge / 2D Corner
428  if (pp.contains("type.xhiylo")) { pp_queryarr("type.xhiylo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XHI_YLO][i] = bcmap[str[i]]; } // 3D Edge / 2D Corner
429  if (pp.contains("type.xhiyhi")) { pp_queryarr("type.xhiyhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XHI_YHI][i] = bcmap[str[i]]; } // 3D Edge / 2D Corner
430 
431  if (pp.contains("type.xlo")) { pp_queryarr("type.xlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XLO][i] = bcmap[str[i]]; } // 3D Face 2D Edge
432  if (pp.contains("type.xhi")) { pp_queryarr("type.xhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::XHI][i] = bcmap[str[i]]; } // 3D Face 2D Edge
433  if (pp.contains("type.ylo")) { pp_queryarr("type.ylo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::YLO][i] = bcmap[str[i]]; } // 3D Face 2D Edge
434  if (pp.contains("type.yhi")) { pp_queryarr("type.yhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::YHI][i] = bcmap[str[i]]; } // 3D Face 2D Edge
435  #if AMREX_SPACEDIM==3
436  if (pp.contains("type.zlo")) { pp_queryarr("type.zlo",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::ZLO][i] = bcmap[str[i]]; } // 3D Face
437  if (pp.contains("type.zhi")) { pp_queryarr("type.zhi",str); for (int i = 0; i < AMREX_SPACEDIM; i++) value.m_bc_type[Face::ZHI][i] = bcmap[str[i]]; } // 3D Face
438  #endif
439 
440  // VALS
441 
442  for (int face = 0; face != Face::INT; face++)
443  {
444  std::vector<std::string> val(AMREX_SPACEDIM,"0.0");
445  val.resize(AMREX_SPACEDIM,"0.0");
446  std::string querystr = std::string("val.") + std::string(facestr[face]);
447  pp_queryarr(querystr.c_str(),val);
448  for (int i = 0 ; i < AMREX_SPACEDIM; i++)
449  {
450  value.m_bc_func_parser[face][i] = amrex::Parser(val[i].c_str());
451  value.m_bc_func_parser[face][i].registerVariables({"x","y","z","t"});
452  value.m_bc_func[face][i] = value.m_bc_func_parser[face][i].compile<4>();
453  }
454  }
455 
456 #if AMREX_SPACEDIM==2
457  // We may wish to use an input file that has 3D BC inputs
458  // This will prevent the parser from complaining that there are unused inputs.
459  std::vector<std::string> ignore_in_2d = {
460  "zlo","zhi",
461  "zhixlo","zloxlo","zhixhi","zloxhi","ylozlo","ylozhi","yhizlo","yhizhi",
462  "xloylozlo","xloylozhi","xloyhizlo","xloyhizhi","xhiylozlo","xhiylozhi","xhiyhizlo","xhiyhizhi"};
463  for (unsigned int n = 0; n < ignore_in_2d.size(); n++)
464  {
465  std::string querystr = std::string("val.") + ignore_in_2d[n];
466  pp.ignore(querystr);
467  querystr = std::string("type.") + ignore_in_2d[n];
468  pp.ignore(querystr);
469  }
470 #endif
471 
472 
473  }
474  Expression(IO::ParmParse &pp, std::string name): Expression()
475  {pp_queryclass(name,*this);}
476 };
477 }
478 }
479 }
480 #endif
Model::Solid::gradu
@ gradu
Definition: Solid.H:26
Set::Position
AMREX_FORCE_INLINE Vector Position(const int &i, const int &j, const int &k, const amrex::Geometry &geom, const amrex::IndexType &ixType)
Definition: Base.H:121
BC::Operator::Elastic::Expression::name
static constexpr const char * name
Definition: Expression.H:20
Linear.H
Operator
Documentation for operator namespace.
Definition: Diagonal.cpp:13
BC::Operator::Elastic::Expression
Definition: Expression.H:17
BC::Operator::Elastic::Expression::Init
virtual void Init(amrex::FabArray< amrex::BaseFab< Set::Vector >> *a_rhs, const amrex::Geometry &a_geom, bool a_homogeneous=false) const override
Definition: Expression.H:54
SQRT2INV
#define SQRT2INV
Definition: Elastic.H:78
BC::Operator::Elastic::Expression::m_bc_func_parser
std::array< std::array< amrex::Parser, AMREX_SPACEDIM >, m_nfaces > m_bc_func_parser
Definition: Expression.H:380
BC::Operator::Elastic::Elastic::m_time
Set::Scalar m_time
Definition: Elastic.H:93
BC::Operator::Elastic::Expression::m_bc_type
std::array< std::array< Type, AMREX_SPACEDIM >, m_nfaces > m_bc_type
Definition: Expression.H:378
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition: Base.H:20
ParmParse.H
Set::None
@ None
Definition: Base.H:197
BC::Operator::Elastic::Elastic::AMREX_D_DECL
@ AMREX_D_DECL
Definition: Elastic.H:44
BC::Operator::Elastic::Elastic::Periodic
@ Periodic
Definition: Elastic.H:24
SQRT3INV
#define SQRT3INV
Definition: Elastic.H:77
pp_queryclass
#define pp_queryclass(...)
Definition: ParmParse.H:106
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition: Base.H:23
Set::Face
@ Face
Definition: Set.H:32
IO::ParmParse::contains
bool contains(std::string name)
Definition: ParmParse.H:153
BC::Operator::Elastic::Expression::Expression
Expression(IO::ParmParse &pp, std::string name)
Definition: Expression.H:474
BC::Operator::Elastic::Elastic
Definition: Elastic.H:18
BC::Operator::Elastic::Expression::getType
virtual std::array< Type, AMREX_SPACEDIM > getType(const int &i, const int &j, [[maybe_unused]] const int &k, const amrex::Box &domain) override
Definition: Expression.H:231
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:166
BC::Operator::Elastic::Expression::~Expression
~Expression()
Definition: Expression.H:48
BC
Collection of boundary condition (BC) objects.
Definition: BC.cpp:4
Elastic.H
IO::ParmParse
Definition: ParmParse.H:112
BC::Operator::Elastic::Elastic::Init
virtual void Init(amrex::MultiFab *a_rhs, const amrex::Geometry &a_geom, bool a_homogeneous=false) const =0
BC::Operator::Elastic::Expression::Init
virtual void Init(amrex::MultiFab *a_rhs, const amrex::Geometry &a_geom, bool a_homogeneous=false) const override
Definition: Expression.H:142
IO::ParmParse::ignore
void ignore(std::string name)
Definition: ParmParse.H:134
BC::Operator::Elastic::Expression::Parse
static void Parse(Expression &value, IO::ParmParse &pp)
Definition: Expression.H:391
BC::Operator::Elastic::Expression::Expression
Expression()
Definition: Expression.H:39
BC::Operator::Elastic::Expression::operator()
AMREX_FORCE_INLINE Set::Vector operator()(const Set::Vector &u, const Set::Matrix &gradu, const Set::Matrix &sigma, const int &i, const int &j, const int &k, const amrex::Box &domain) override
Definition: Expression.H:288
INFO
#define INFO
Definition: Util.H:20
BC::Operator::Elastic::Expression::m_bc_func
std::array< std::array< amrex::ParserExecutor< 4 >, AMREX_SPACEDIM >, m_nfaces > m_bc_func
Definition: Expression.H:381
BC::Operator::Elastic::Expression::set
AMREX_FORCE_INLINE Set::Vector set(std::array< Type, AMREX_SPACEDIM > &bc_type, const Set::Vector &u, const Set::Matrix &gradu, const Set::Matrix &sigma, Set::Vector n) const
Definition: Expression.H:352
pp_queryarr
#define pp_queryarr(...)
Definition: ParmParse.H:103
Generated on Sat Feb 8 2025 00:07:39 for Alamo by   doxygen 1.8.17