1 #ifndef BC_OPERATOR_ELASTIC_EXPRESSION_H
2 #define BC_OPERATOR_ELASTIC_EXPRESSION_H
9 #include "AMReX_Parser.H"
22 static const constexpr
char *
const facestr[] = {
23 "xlo",
"ylo",
"xhi",
"yhi",
"xloylo",
"xloyhi",
"xhiylo",
"xhiyhi",
"int"
25 #elif AMREX_SPACEDIM==3
26 static const constexpr
char *
const facestr[] = {
27 "xlo",
"ylo",
"zlo",
"xhi",
"yhi",
"zhi",
28 "ylozlo",
"ylozhi",
"yhizlo",
"yhizhi",
29 "zloxlo",
"zloxhi",
"zhixlo",
"zhixhi",
30 "xloylo",
"xloyhi",
"xhiylo",
"xhiyhi",
31 "xloylozlo",
"xloylozhi",
"xloyhizlo",
"xloyhizhi",
32 "xhiylozlo",
"xhiylozhi",
"xhiyhizlo",
"xhiyhizhi",
40 for (
int face = 0; face < m_nfaces; face++)
41 for (
int direction = 0; direction < AMREX_SPACEDIM; direction++)
43 m_bc_type[face][direction] = Type::Displacement;
52 Init(amrex::FabArray<amrex::BaseFab<Set::Vector>> *a_rhs,
53 const amrex::Geometry &a_geom,
54 bool a_homogeneous =
false)
const override
56 amrex::Box domain(a_geom.Domain());
57 domain.convert(amrex::IntVect::TheNodeVector());
58 const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
59 for (amrex::MFIter mfi(*a_rhs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
61 amrex::Box bx = mfi.tilebox();
64 amrex::Array4<Set::Vector>
const& rhs = a_rhs->array(mfi);
65 amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(
int i,
int j,
int k) {
67 for (
int dir = 0; dir < AMREX_SPACEDIM; dir++)
69 Face face = Face::INT;
71 #if AMREX_SPACEDIM == 2
73 if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
74 else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
75 else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
76 else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
78 else if (i==lo.x) face = Face::XLO;
79 else if (i==hi.x) face = Face::XHI;
80 else if (j==lo.y) face = Face::YLO;
81 else if (j==hi.y) face = Face::YHI;
83 #elif AMREX_SPACEDIM == 3
85 if (i==lo.x && j==lo.y && k==lo.z) face = Face::XLO_YLO_ZLO;
86 else if (i==lo.x && j==lo.y && k==hi.z) face = Face::XLO_YLO_ZHI;
87 else if (i==lo.x && j==hi.y && k==lo.z) face = Face::XLO_YHI_ZLO;
88 else if (i==lo.x && j==hi.y && k==hi.z) face = Face::XLO_YHI_ZHI;
89 else if (i==hi.x && j==lo.y && k==lo.z) face = Face::XHI_YLO_ZLO;
90 else if (i==hi.x && j==lo.y && k==hi.z) face = Face::XHI_YLO_ZHI;
91 else if (i==hi.x && j==hi.y && k==lo.z) face = Face::XHI_YHI_ZLO;
92 else if (i==hi.x && j==hi.y && k==hi.z) face = Face::XHI_YHI_ZHI;
94 else if (j==lo.y && k==lo.z) face = Face::YLO_ZLO;
95 else if (j==lo.y && k==hi.z) face = Face::YLO_ZHI;
96 else if (j==hi.y && k==lo.z) face = Face::YHI_ZLO;
97 else if (j==hi.y && k==hi.z) face = Face::YHI_ZHI;
98 else if (k==lo.z && i==lo.x) face = Face::ZLO_XLO;
99 else if (k==lo.z && i==hi.x) face = Face::ZLO_XHI;
100 else if (k==hi.z && i==lo.x) face = Face::ZHI_XLO;
101 else if (k==hi.z && i==hi.x) face = Face::ZHI_XHI;
102 else if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
103 else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
104 else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
105 else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
107 else if (i==lo.x) face = Face::XLO;
108 else if (i==hi.x) face = Face::XHI;
109 else if (j==lo.y) face = Face::YLO;
110 else if (j==hi.y) face = Face::YHI;
111 else if (k==lo.z) face = Face::ZLO;
112 else if (k==hi.z) face = Face::ZHI;
116 amrex::IndexType type = a_rhs->ixType();
117 if (!(face == Face::INT))
119 if (a_homogeneous &&
m_bc_type[face][dir] == Type::Displacement)
120 rhs(i,j,k)(dir) = 0.0;
124 #if AMREX_SPACEDIM==1
126 #elif AMREX_SPACEDIM==2
128 #elif AMREX_SPACEDIM==3
141 const amrex::Geometry &a_geom,
142 bool a_homogeneous =
false)
const override
144 amrex::Box domain(a_geom.Domain());
145 domain.convert(amrex::IntVect::TheNodeVector());
146 const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
147 for (amrex::MFIter mfi(*a_rhs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
149 amrex::Box bx = mfi.tilebox();
152 amrex::Array4<amrex::Real>
const& rhs = a_rhs->array(mfi);
153 amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(
int i,
int j,
int k) {
155 for (
int dir = 0; dir < AMREX_SPACEDIM; dir++)
157 Face face = Face::INT;
159 #if AMREX_SPACEDIM == 2
161 if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
162 else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
163 else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
164 else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
166 else if (i==lo.x) face = Face::XLO;
167 else if (i==hi.x) face = Face::XHI;
168 else if (j==lo.y) face = Face::YLO;
169 else if (j==hi.y) face = Face::YHI;
171 #elif AMREX_SPACEDIM == 3
173 if (i==lo.x && j==lo.y && k==lo.z) face = Face::XLO_YLO_ZLO;
174 else if (i==lo.x && j==lo.y && k==hi.z) face = Face::XLO_YLO_ZHI;
175 else if (i==lo.x && j==hi.y && k==lo.z) face = Face::XLO_YHI_ZLO;
176 else if (i==lo.x && j==hi.y && k==hi.z) face = Face::XLO_YHI_ZHI;
177 else if (i==hi.x && j==lo.y && k==lo.z) face = Face::XHI_YLO_ZLO;
178 else if (i==hi.x && j==lo.y && k==hi.z) face = Face::XHI_YLO_ZHI;
179 else if (i==hi.x && j==hi.y && k==lo.z) face = Face::XHI_YHI_ZLO;
180 else if (i==hi.x && j==hi.y && k==hi.z) face = Face::XHI_YHI_ZHI;
182 else if (j==lo.y && k==lo.z) face = Face::YLO_ZLO;
183 else if (j==lo.y && k==hi.z) face = Face::YLO_ZHI;
184 else if (j==hi.y && k==lo.z) face = Face::YHI_ZLO;
185 else if (j==hi.y && k==hi.z) face = Face::YHI_ZHI;
186 else if (k==lo.z && i==lo.x) face = Face::ZLO_XLO;
187 else if (k==lo.z && i==hi.x) face = Face::ZLO_XHI;
188 else if (k==hi.z && i==lo.x) face = Face::ZHI_XLO;
189 else if (k==hi.z && i==hi.x) face = Face::ZHI_XHI;
190 else if (i==lo.x && j==lo.y) face = Face::XLO_YLO;
191 else if (i==lo.x && j==hi.y) face = Face::XLO_YHI;
192 else if (i==hi.x && j==lo.y) face = Face::XHI_YLO;
193 else if (i==hi.x && j==hi.y) face = Face::XHI_YHI;
195 else if (i==lo.x) face = Face::XLO;
196 else if (i==hi.x) face = Face::XHI;
197 else if (j==lo.y) face = Face::YLO;
198 else if (j==hi.y) face = Face::YHI;
199 else if (k==lo.z) face = Face::ZLO;
200 else if (k==hi.z) face = Face::ZHI;
204 amrex::IndexType type = a_rhs->ixType();
205 if (!(face == Face::INT))
207 if (a_homogeneous &&
m_bc_type[face][dir] == Type::Displacement)
208 rhs(i,j,k,dir) = 0.0;
212 #if AMREX_SPACEDIM==1
214 #elif AMREX_SPACEDIM==2
216 #elif AMREX_SPACEDIM==3
230 const int &i,
const int &j, [[maybe_unused]]
const int &k,
231 const amrex::Box &domain)
override
233 (void)i, (
void)j, (void)k;
235 const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
238 #if AMREX_SPACEDIM == 2
239 if (m[0] == lo.x && m[1] == lo.y)
return m_bc_type[Face::XLO_YLO];
240 if (m[0] == lo.x && m[1] == hi.y)
return m_bc_type[Face::XLO_YHI];
241 if (m[0] == hi.x && m[1] == lo.y)
return m_bc_type[Face::XHI_YLO];
242 if (m[0] == hi.x && m[1] == hi.y)
return m_bc_type[Face::XHI_YHI];
243 if (m[0] == lo.x)
return m_bc_type[Face::XLO];
244 if (m[0] == hi.x)
return m_bc_type[Face::XHI];
245 if (m[1] == lo.y)
return m_bc_type[Face::YLO];
246 if (m[1] == hi.y)
return m_bc_type[Face::YHI];
248 #elif AMREX_SPACEDIM == 3
250 if (m[0] == lo.x && m[1] == lo.y && m[2] == lo.z)
return m_bc_type[Face::XLO_YLO_ZLO];
251 if (m[0] == lo.x && m[1] == lo.y && m[2] == hi.z)
return m_bc_type[Face::XLO_YLO_ZHI];
252 if (m[0] == lo.x && m[1] == hi.y && m[2] == lo.z)
return m_bc_type[Face::XLO_YHI_ZLO];
253 if (m[0] == lo.x && m[1] == hi.y && m[2] == hi.z)
return m_bc_type[Face::XLO_YHI_ZHI];
254 if (m[0] == hi.x && m[1] == lo.y && m[2] == lo.z)
return m_bc_type[Face::XHI_YLO_ZLO];
255 if (m[0] == hi.x && m[1] == lo.y && m[2] == hi.z)
return m_bc_type[Face::XHI_YLO_ZHI];
256 if (m[0] == hi.x && m[1] == hi.y && m[2] == lo.z)
return m_bc_type[Face::XHI_YHI_ZLO];
257 if (m[0] == hi.x && m[1] == hi.y && m[2] == hi.z)
return m_bc_type[Face::XHI_YHI_ZHI];
259 if (m[1] == lo.y && m[2] == lo.z)
return m_bc_type[Face::YLO_ZLO];
260 if (m[1] == lo.y && m[2] == hi.z)
return m_bc_type[Face::YLO_ZHI];
261 if (m[1] == hi.y && m[2] == lo.z)
return m_bc_type[Face::YHI_ZLO];
262 if (m[1] == hi.y && m[2] == hi.z)
return m_bc_type[Face::YHI_ZHI];
263 if (m[2] == lo.z && m[0] == lo.x)
return m_bc_type[Face::ZLO_XLO];
264 if (m[2] == lo.z && m[0] == hi.x)
return m_bc_type[Face::ZLO_XHI];
265 if (m[2] == hi.z && m[0] == lo.x)
return m_bc_type[Face::ZHI_XLO];
266 if (m[2] == hi.z && m[0] == hi.x)
return m_bc_type[Face::ZHI_XHI];
267 if (m[0] == lo.x && m[1] == lo.y)
return m_bc_type[Face::XLO_YLO];
268 if (m[0] == lo.x && m[1] == hi.y)
return m_bc_type[Face::XLO_YHI];
269 if (m[0] == hi.x && m[1] == lo.y)
return m_bc_type[Face::XHI_YLO];
270 if (m[0] == hi.x && m[1] == hi.y)
return m_bc_type[Face::XHI_YHI];
272 if (m[0] == lo.x)
return m_bc_type[Face::XLO];
273 if (m[0] == hi.x)
return m_bc_type[Face::XHI];
274 if (m[1] == lo.y)
return m_bc_type[Face::YLO];
275 if (m[1] == hi.y)
return m_bc_type[Face::YHI];
276 if (m[2] == lo.z)
return m_bc_type[Face::ZLO];
277 if (m[2] == hi.z)
return m_bc_type[Face::ZHI];
289 const int &i,
const int &j,
const int &k,
290 const amrex::Box &domain)
override
292 (void)i; (void)j; (void)k;
295 const amrex::Dim3 lo= amrex::lbound(domain), hi = amrex::ubound(domain);
300 #if AMREX_SPACEDIM == 2
312 #elif AMREX_SPACEDIM == 3
346 return Set::Vector::Zero();
354 for (
int i = 0; i < AMREX_SPACEDIM; i++)
356 if (bc_type[i] == Type::Displacement)
358 else if (bc_type[i] == Type::Traction)
360 else if (bc_type[i] == Type::Neumann)
370 #if AMREX_SPACEDIM==2
371 static const int m_nfaces = 8;
372 #elif AMREX_SPACEDIM==3
373 static const int m_nfaces = 26;
376 std::array<std::array<Type, AMREX_SPACEDIM>, m_nfaces>
m_bc_type;
379 std::array<std::array<amrex::ParserExecutor<4>,AMREX_SPACEDIM>, m_nfaces>
m_bc_func;
391 std::map<std::string, Type> bcmap;
392 bcmap[
"displacement"] = Type::Displacement;
393 bcmap[
"disp"] = Type::Displacement;
394 bcmap[
"neumann"] = Type::Neumann;
395 bcmap[
"traction"] = Type::Traction;
396 bcmap[
"trac"] = Type::Traction;
397 bcmap[
"periodic"] = Type::Periodic;
399 std::vector<std::string> str;
403 #if AMREX_SPACEDIM==3
404 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]]; }
405 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]]; }
406 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]]; }
407 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]]; }
408 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]]; }
409 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]]; }
410 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]]; }
411 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]]; }
414 #if AMREX_SPACEDIM==3
415 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]]; }
416 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]]; }
417 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]]; }
418 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]]; }
419 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]]; }
420 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]]; }
421 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]]; }
422 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]]; }
424 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]]; }
425 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]]; }
426 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]]; }
427 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]]; }
429 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]]; }
430 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]]; }
431 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]]; }
432 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]]; }
433 #if AMREX_SPACEDIM==3
434 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]]; }
435 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]]; }
440 for (
int face = 0; face != Face::INT; face++)
442 std::vector<std::string> val(AMREX_SPACEDIM,
"0.0");
443 val.resize(AMREX_SPACEDIM,
"0.0");
444 std::string querystr = std::string(
"val.") + std::string(facestr[face]);
446 for (
int i = 0 ; i < AMREX_SPACEDIM; i++)
454 #if AMREX_SPACEDIM==2
457 std::vector<std::string> ignore_in_2d = {
459 "zhixlo",
"zloxlo",
"zhixhi",
"zloxhi",
"ylozlo",
"ylozhi",
"yhizlo",
"yhizhi",
460 "xloylozlo",
"xloylozhi",
"xloyhizlo",
"xloyhizhi",
"xhiylozlo",
"xhiylozhi",
"xhiyhizlo",
"xhiyhizhi"};
461 for (
unsigned int n = 0; n < ignore_in_2d.size(); n++)
463 std::string querystr = std::string(
"val.") + ignore_in_2d[n];
465 querystr = std::string(
"type.") + ignore_in_2d[n];