Line data Source code
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 1 : Expression()
40 1 : {
41 : // By default, all boundary conditions are displacement
42 9 : for (int face = 0; face < m_nfaces; face++)
43 24 : for (int direction = 0; direction < AMREX_SPACEDIM; direction++)
44 : {
45 16 : m_bc_type[face][direction] = Type::Displacement;
46 : }
47 1 : };
48 0 : ~Expression() {};
49 :
50 :
51 : using Elastic::Init;
52 :
53 : virtual void
54 2 : Init(amrex::FabArray<amrex::BaseFab<Set::Vector>> *a_rhs,
55 : const amrex::Geometry &a_geom,
56 : bool a_homogeneous = false) const override
57 : {
58 2 : amrex::Box domain(a_geom.Domain());
59 4 : domain.convert(amrex::IntVect::TheNodeVector());
60 2 : const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
61 4 : for (amrex::MFIter mfi(*a_rhs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
62 : {
63 2 : amrex::Box bx = mfi.tilebox();
64 2 : bx.grow(2);
65 2 : bx = bx & domain;
66 2 : amrex::Array4<Set::Vector> const& rhs = a_rhs->array(mfi);
67 2 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
68 :
69 1734 : for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
70 : {
71 1156 : Face face = Face::INT;
72 :
73 : #if AMREX_SPACEDIM == 2
74 :
75 1156 : if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
76 1152 : else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
77 1148 : else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
78 1144 : else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
79 :
80 1140 : else if (i==lo.x) face = Face::XLO;
81 1080 : else if (i==hi.x) face = Face::XHI;
82 1020 : else if (j==lo.y) face = Face::YLO;
83 960 : else if (j==hi.y) face = Face::YHI;
84 :
85 : #elif AMREX_SPACEDIM == 3
86 :
87 0 : if (i==lo.x && j==lo.y && k==lo.z) face = Face::XLO_YLO_ZLO;
88 0 : else if (i==lo.x && j==lo.y && k==hi.z) face = Face::XLO_YLO_ZHI;
89 0 : else if (i==lo.x && j==hi.y && k==lo.z) face = Face::XLO_YHI_ZLO;
90 0 : else if (i==lo.x && j==hi.y && k==hi.z) face = Face::XLO_YHI_ZHI;
91 0 : else if (i==hi.x && j==lo.y && k==lo.z) face = Face::XHI_YLO_ZLO;
92 0 : else if (i==hi.x && j==lo.y && k==hi.z) face = Face::XHI_YLO_ZHI;
93 0 : else if (i==hi.x && j==hi.y && k==lo.z) face = Face::XHI_YHI_ZLO;
94 0 : else if (i==hi.x && j==hi.y && k==hi.z) face = Face::XHI_YHI_ZHI;
95 :
96 0 : else if (j==lo.y && k==lo.z) face = Face::YLO_ZLO;
97 0 : else if (j==lo.y && k==hi.z) face = Face::YLO_ZHI;
98 0 : else if (j==hi.y && k==lo.z) face = Face::YHI_ZLO;
99 0 : else if (j==hi.y && k==hi.z) face = Face::YHI_ZHI;
100 0 : else if (k==lo.z && i==lo.x) face = Face::ZLO_XLO;
101 0 : else if (k==lo.z && i==hi.x) face = Face::ZLO_XHI;
102 0 : else if (k==hi.z && i==lo.x) face = Face::ZHI_XLO;
103 0 : else if (k==hi.z && i==hi.x) face = Face::ZHI_XHI;
104 0 : else if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
105 0 : else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
106 0 : else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
107 0 : else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
108 :
109 0 : else if (i==lo.x) face = Face::XLO;
110 0 : else if (i==hi.x) face = Face::XHI;
111 0 : else if (j==lo.y) face = Face::YLO;
112 0 : else if (j==hi.y) face = Face::YHI;
113 0 : else if (k==lo.z) face = Face::ZLO;
114 0 : else if (k==hi.z) face = Face::ZHI;
115 :
116 : #endif
117 :
118 1156 : amrex::IndexType type = a_rhs->ixType();
119 1156 : if (!(face == Face::INT))
120 : {
121 256 : if (a_homogeneous && m_bc_type[face][dir] == Type::Displacement)
122 0 : rhs(i,j,k)(dir) = 0.0;
123 : else
124 : {
125 256 : 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 768 : rhs(i,j,k)(dir) = m_bc_func[face][dir](x(0),x(1),0.0,m_time);
130 : #elif AMREX_SPACEDIM==3
131 0 : rhs(i,j,k)(dir) = m_bc_func[face][dir](x(0),x(1),x(2),m_time);
132 : #endif
133 : }
134 : }
135 : }
136 578 : });
137 : }
138 2 : }
139 :
140 :
141 : virtual void
142 0 : Init(amrex::MultiFab * a_rhs,
143 : const amrex::Geometry &a_geom,
144 : bool a_homogeneous = false) const override
145 : {
146 0 : amrex::Box domain(a_geom.Domain());
147 0 : domain.convert(amrex::IntVect::TheNodeVector());
148 0 : const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
149 0 : for (amrex::MFIter mfi(*a_rhs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
150 : {
151 0 : amrex::Box bx = mfi.tilebox();
152 0 : bx.grow(2);
153 0 : bx = bx & domain;
154 0 : amrex::Array4<amrex::Real> const& rhs = a_rhs->array(mfi);
155 0 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
156 :
157 0 : for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
158 : {
159 0 : Face face = Face::INT;
160 :
161 : #if AMREX_SPACEDIM == 2
162 :
163 0 : if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
164 0 : else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
165 0 : else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
166 0 : else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
167 :
168 0 : else if (i==lo.x) face = Face::XLO;
169 0 : else if (i==hi.x) face = Face::XHI;
170 0 : else if (j==lo.y) face = Face::YLO;
171 0 : else if (j==hi.y) face = Face::YHI;
172 :
173 : #elif AMREX_SPACEDIM == 3
174 :
175 0 : if (i==lo.x && j==lo.y && k==lo.z) face = Face::XLO_YLO_ZLO;
176 0 : else if (i==lo.x && j==lo.y && k==hi.z) face = Face::XLO_YLO_ZHI;
177 0 : else if (i==lo.x && j==hi.y && k==lo.z) face = Face::XLO_YHI_ZLO;
178 0 : else if (i==lo.x && j==hi.y && k==hi.z) face = Face::XLO_YHI_ZHI;
179 0 : else if (i==hi.x && j==lo.y && k==lo.z) face = Face::XHI_YLO_ZLO;
180 0 : else if (i==hi.x && j==lo.y && k==hi.z) face = Face::XHI_YLO_ZHI;
181 0 : else if (i==hi.x && j==hi.y && k==lo.z) face = Face::XHI_YHI_ZLO;
182 0 : else if (i==hi.x && j==hi.y && k==hi.z) face = Face::XHI_YHI_ZHI;
183 :
184 0 : else if (j==lo.y && k==lo.z) face = Face::YLO_ZLO;
185 0 : else if (j==lo.y && k==hi.z) face = Face::YLO_ZHI;
186 0 : else if (j==hi.y && k==lo.z) face = Face::YHI_ZLO;
187 0 : else if (j==hi.y && k==hi.z) face = Face::YHI_ZHI;
188 0 : else if (k==lo.z && i==lo.x) face = Face::ZLO_XLO;
189 0 : else if (k==lo.z && i==hi.x) face = Face::ZLO_XHI;
190 0 : else if (k==hi.z && i==lo.x) face = Face::ZHI_XLO;
191 0 : else if (k==hi.z && i==hi.x) face = Face::ZHI_XHI;
192 0 : else if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
193 0 : else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
194 0 : else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
195 0 : else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
196 :
197 0 : else if (i==lo.x) face = Face::XLO;
198 0 : else if (i==hi.x) face = Face::XHI;
199 0 : else if (j==lo.y) face = Face::YLO;
200 0 : else if (j==hi.y) face = Face::YHI;
201 0 : else if (k==lo.z) face = Face::ZLO;
202 0 : else if (k==hi.z) face = Face::ZHI;
203 :
204 : #endif
205 :
206 0 : amrex::IndexType type = a_rhs->ixType();
207 0 : if (!(face == Face::INT))
208 : {
209 0 : if (a_homogeneous && m_bc_type[face][dir] == Type::Displacement)
210 0 : rhs(i,j,k,dir) = 0.0;
211 : else
212 : {
213 0 : 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 0 : rhs(i,j,k,dir) = (m_bc_func[face][dir])(x(0),x(1),0.0,m_time);
218 : #elif AMREX_SPACEDIM==3
219 0 : rhs(i,j,k,dir) = (m_bc_func[face][dir])(x(0),x(1),x(2),m_time);
220 : #endif
221 : }
222 : }
223 : }
224 0 : });
225 : }
226 0 : }
227 :
228 :
229 :
230 : virtual
231 0 : 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 0 : 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 0 : if (m[0] == lo.x && m[1] == lo.y) return m_bc_type[Face::XLO_YLO];
242 0 : if (m[0] == lo.x && m[1] == hi.y) return m_bc_type[Face::XLO_YHI];
243 0 : if (m[0] == hi.x && m[1] == lo.y) return m_bc_type[Face::XHI_YLO];
244 0 : if (m[0] == hi.x && m[1] == hi.y) return m_bc_type[Face::XHI_YHI];
245 0 : if (m[0] == lo.x) return m_bc_type[Face::XLO];
246 0 : if (m[0] == hi.x) return m_bc_type[Face::XHI];
247 0 : if (m[1] == lo.y) return m_bc_type[Face::YLO];
248 0 : if (m[1] == hi.y) return m_bc_type[Face::YHI];
249 :
250 : #elif AMREX_SPACEDIM == 3
251 :
252 0 : if (m[0] == lo.x && m[1] == lo.y && m[2] == lo.z) return m_bc_type[Face::XLO_YLO_ZLO];
253 0 : if (m[0] == lo.x && m[1] == lo.y && m[2] == hi.z) return m_bc_type[Face::XLO_YLO_ZHI];
254 0 : if (m[0] == lo.x && m[1] == hi.y && m[2] == lo.z) return m_bc_type[Face::XLO_YHI_ZLO];
255 0 : if (m[0] == lo.x && m[1] == hi.y && m[2] == hi.z) return m_bc_type[Face::XLO_YHI_ZHI];
256 0 : if (m[0] == hi.x && m[1] == lo.y && m[2] == lo.z) return m_bc_type[Face::XHI_YLO_ZLO];
257 0 : if (m[0] == hi.x && m[1] == lo.y && m[2] == hi.z) return m_bc_type[Face::XHI_YLO_ZHI];
258 0 : if (m[0] == hi.x && m[1] == hi.y && m[2] == lo.z) return m_bc_type[Face::XHI_YHI_ZLO];
259 0 : if (m[0] == hi.x && m[1] == hi.y && m[2] == hi.z) return m_bc_type[Face::XHI_YHI_ZHI];
260 :
261 0 : if (m[1] == lo.y && m[2] == lo.z) return m_bc_type[Face::YLO_ZLO];
262 0 : if (m[1] == lo.y && m[2] == hi.z) return m_bc_type[Face::YLO_ZHI];
263 0 : if (m[1] == hi.y && m[2] == lo.z) return m_bc_type[Face::YHI_ZLO];
264 0 : if (m[1] == hi.y && m[2] == hi.z) return m_bc_type[Face::YHI_ZHI];
265 0 : if (m[2] == lo.z && m[0] == lo.x) return m_bc_type[Face::ZLO_XLO];
266 0 : if (m[2] == lo.z && m[0] == hi.x) return m_bc_type[Face::ZLO_XHI];
267 0 : if (m[2] == hi.z && m[0] == lo.x) return m_bc_type[Face::ZHI_XLO];
268 0 : if (m[2] == hi.z && m[0] == hi.x) return m_bc_type[Face::ZHI_XHI];
269 0 : if (m[0] == lo.x && m[1] == lo.y) return m_bc_type[Face::XLO_YLO];
270 0 : if (m[0] == lo.x && m[1] == hi.y) return m_bc_type[Face::XLO_YHI];
271 0 : if (m[0] == hi.x && m[1] == lo.y) return m_bc_type[Face::XHI_YLO];
272 0 : if (m[0] == hi.x && m[1] == hi.y) return m_bc_type[Face::XHI_YHI];
273 :
274 0 : if (m[0] == lo.x) return m_bc_type[Face::XLO];
275 0 : if (m[0] == hi.x) return m_bc_type[Face::XHI];
276 0 : if (m[1] == lo.y) return m_bc_type[Face::YLO];
277 0 : if (m[1] == hi.y) return m_bc_type[Face::YHI];
278 0 : if (m[2] == lo.z) return m_bc_type[Face::ZLO];
279 0 : if (m[2] == hi.z) return m_bc_type[Face::ZHI];
280 :
281 : #endif
282 :
283 0 : return {AMREX_D_DECL(Type::None,Type::None,Type::None)};
284 : }
285 :
286 :
287 : AMREX_FORCE_INLINE
288 58256 : 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 58256 : amrex::IntVect m(AMREX_D_DECL(i,j,k));
300 :
301 : // Corners
302 : #if AMREX_SPACEDIM == 2
303 :
304 77272 : 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 72820 : 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 72820 : 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 68368 : 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 61690 : if (m[0] == lo.x) return set(m_bc_type[Face::XLO], u, gradu, sigma, Set::Vector(-1, 0));
310 49352 : if (m[0] == hi.x) return set(m_bc_type[Face::XHI], u, gradu, sigma, Set::Vector(+1, 0));
311 37014 : if (m[1] == lo.y) return set(m_bc_type[Face::YLO], u, gradu, sigma, Set::Vector( 0,-1));
312 24676 : 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 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));
317 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));
318 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));
319 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));
320 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));
321 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));
322 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));
323 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));
324 :
325 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));
326 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));
327 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));
328 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));
329 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));
330 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));
331 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));
332 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));
333 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 ));
334 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 ));
335 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 ));
336 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 ));
337 :
338 0 : if (m[0] == lo.x) return set(m_bc_type[Face::XLO], u, gradu, sigma, Set::Vector(-1, 0, 0));
339 0 : if (m[0] == hi.x) return set(m_bc_type[Face::XHI], u, gradu, sigma, Set::Vector(+1, 0, 0));
340 0 : if (m[1] == lo.y) return set(m_bc_type[Face::YLO], u, gradu, sigma, Set::Vector( 0,-1, 0));
341 0 : if (m[1] == hi.y) return set(m_bc_type[Face::YHI], u, gradu, sigma, Set::Vector( 0,+1, 0));
342 0 : if (m[2] == lo.z) return set(m_bc_type[Face::ZLO], u, gradu, sigma, Set::Vector( 0, 0,-1));
343 0 : 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 0 : Util::Abort(INFO,"Boundary condition error");
348 0 : 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 58256 : Set::Vector f = Set::Vector::Zero();
356 174768 : for (int i = 0; i < AMREX_SPACEDIM; i++)
357 : {
358 116512 : if (bc_type[i] == Type::Displacement)
359 91836 : f(i) = u(i);
360 24676 : else if (bc_type[i] == Type::Traction)
361 24676 : f(i) = (sigma*n)(i);
362 0 : else if (bc_type[i] == Type::Neumann)
363 0 : f(i) = (gradu*n)(i);
364 0 : else if (bc_type[i] == Periodic)
365 0 : continue;
366 : }
367 58256 : 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 1 : static void Parse(Expression & value, IO::ParmParse & pp)
392 : {
393 2 : std::map<std::string, Type> bcmap;
394 1 : bcmap["displacement"] = Type::Displacement;
395 1 : bcmap["disp"] = Type::Displacement;
396 1 : bcmap["neumann"] = Type::Neumann;
397 1 : bcmap["traction"] = Type::Traction;
398 1 : bcmap["trac"] = Type::Traction;
399 1 : bcmap["periodic"] = Type::Periodic;
400 :
401 2 : std::vector<std::string> str;
402 :
403 : // TYPES
404 :
405 : #if AMREX_SPACEDIM==3
406 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 1 : 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 1 : 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 1 : 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 1 : 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 1 : 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 1 : 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 1 : 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 3 : 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 0 : 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 0 : 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 9 : for (int face = 0; face != Face::INT; face++)
443 : {
444 24 : std::vector<std::string> val(AMREX_SPACEDIM,"0.0");
445 8 : val.resize(AMREX_SPACEDIM,"0.0");
446 24 : std::string querystr = std::string("val.") + std::string(facestr[face]);
447 8 : pp_queryarr(querystr.c_str(),val);
448 24 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
449 : {
450 16 : value.m_bc_func_parser[face][i] = amrex::Parser(val[i].c_str());
451 80 : value.m_bc_func_parser[face][i].registerVariables({"x","y","z","t"});
452 16 : 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 22 : "xloylozlo","xloylozhi","xloyhizlo","xloyhizhi","xhiylozlo","xhiylozhi","xhiyhizlo","xhiyhizhi"};
463 19 : for (unsigned int n = 0; n < ignore_in_2d.size(); n++)
464 : {
465 36 : std::string querystr = std::string("val.") + ignore_in_2d[n];
466 18 : pp.ignore(querystr);
467 18 : querystr = std::string("type.") + ignore_in_2d[n];
468 18 : pp.ignore(querystr);
469 : }
470 : #endif
471 :
472 :
473 1 : }
474 1 : Expression(IO::ParmParse &pp, std::string name): Expression()
475 1 : {pp_queryclass(name,*this);}
476 : };
477 : }
478 : }
479 : }
480 : #endif
|
