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