Alamo
Elastic.cpp
Go to the documentation of this file.
1// TODO: Remove these
2
3#include "Elastic.H"
4#include "Set/Set.H"
5
6#include "Numeric/Stencil.H"
7namespace Operator
8{
9template<int SYM>
10Elastic<SYM>::Elastic(const Vector<Geometry>& a_geom,
11 const Vector<BoxArray>& a_grids,
12 const Vector<DistributionMapping>& a_dmap,
13 const LPInfo& a_info)
14{
15 BL_PROFILE("Operator::Elastic::Elastic()");
16
17 define(a_geom, a_grids, a_dmap, a_info);
18}
19
20template<int SYM>
23
24template<int SYM>
25void
26Elastic<SYM>::define(const Vector<Geometry>& a_geom,
27 const Vector<BoxArray>& a_grids,
28 const Vector<DistributionMapping>& a_dmap,
29 const LPInfo& a_info,
30 const Vector<FabFactory<FArrayBox> const*>& a_factory)
31{
32 BL_PROFILE("Operator::Elastic::define()");
33
34 Operator::define(a_geom, a_grids, a_dmap, a_info, a_factory);
35
36 int model_nghost = 2;
37
38 m_ddw_mf.resize(m_num_amr_levels);
39 m_psi_mf.resize(m_num_amr_levels);
40 for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
41 {
42 m_ddw_mf[amrlev].resize(m_num_mg_levels[amrlev]);
43 m_psi_mf[amrlev].resize(m_num_mg_levels[amrlev]);
44 for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
45 {
46 m_ddw_mf[amrlev][mglev].reset(new MultiTab(amrex::convert(m_grids[amrlev][mglev],
47 amrex::IntVect::TheNodeVector()),
48 m_dmap[amrlev][mglev], 1, model_nghost));
49 m_psi_mf[amrlev][mglev].reset(new MultiFab(m_grids[amrlev][mglev],
50 m_dmap[amrlev][mglev], 1, model_nghost));
51
52 if (!m_psi_set) m_psi_mf[amrlev][mglev]->setVal(1.0);
53 }
54 }
55}
56
57template <int SYM>
58void
59Elastic<SYM>::SetModel(MATRIX4& a_model)
60{
61 for (int amrlev = 0; amrlev < m_num_amr_levels; amrlev++)
62 {
63 amrex::Box domain(m_geom[amrlev][0].Domain());
64 domain.convert(amrex::IntVect::TheNodeVector());
65
66 for (MFIter mfi(*m_ddw_mf[amrlev][0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
67 {
68 Box bx = mfi.grownnodaltilebox();
69
70 amrex::Array4<MATRIX4> const& ddw = (*(m_ddw_mf[amrlev][0])).array(mfi);
71
72 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
73 ddw(i, j, k) = a_model;
74
75#ifdef AMREX_DEBUG
76 if (ddw(i, j, k).contains_nan()) Util::Abort(INFO, "model is nan at (", i, ",", j, ",", k, "), amrlev=", amrlev);
77#endif
78 });
79 }
80
81 (*(m_ddw_mf[amrlev][0])).setMultiGhost(true);
82 (*(m_ddw_mf[amrlev][0])).FillBoundary( Geom(amrlev,0).periodicity());
83 }
84 m_model_set = true;
85}
86
87template <int SYM>
88void
89Elastic<SYM>::SetModel(int amrlev, const amrex::FabArray<amrex::BaseFab<MATRIX4> >& a_model)
90{
91 BL_PROFILE("Operator::Elastic::SetModel()");
92
93 amrex::Box domain(m_geom[amrlev][0].Domain());
94 domain.convert(amrex::IntVect::TheNodeVector());
95
96 if (a_model.boxArray() != m_ddw_mf[amrlev][0]->boxArray()) Util::Abort(INFO, "Inconsistent box arrays\n", "a_model.boxArray()=\n", a_model.boxArray(), "\n but the current box array is \n", m_ddw_mf[amrlev][0]->boxArray());
97 if (a_model.DistributionMap() != m_ddw_mf[amrlev][0]->DistributionMap()) Util::Abort(INFO, "Inconsistent distribution maps");
98 if (a_model.nComp() != m_ddw_mf[amrlev][0]->nComp()) Util::Abort(INFO, "Inconsistent # of components - should be ", m_ddw_mf[amrlev][0]->nComp());
99 if (a_model.nGrow() != m_ddw_mf[amrlev][0]->nGrow()) Util::Abort(INFO, "Inconsistent # of ghost nodes, should be ", m_ddw_mf[amrlev][0]->nGrow());
100
101
102 for (MFIter mfi(a_model, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
103 {
104 Box bx = mfi.grownnodaltilebox();
105
106 amrex::Array4<MATRIX4> const& C = (*(m_ddw_mf[amrlev][0])).array(mfi);
107 amrex::Array4<const MATRIX4> const& a_C = a_model.array(mfi);
108
109 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
110 C(i, j, k) = a_C(i, j, k);
111 });
112 }
113 m_ddw_mf[amrlev][0]->setMultiGhost(true);
114 m_ddw_mf[amrlev][0]->FillBoundaryAndSync(Geom(amrlev,0).periodicity());
115 m_model_set = true;
116}
117
118template <int SYM>
119void
120Elastic<SYM>::SetPsi(int amrlev, const amrex::MultiFab& a_psi_mf)
121{
122 BL_PROFILE("Operator::Elastic::SetPsi()");
123 amrex::Box domain(m_geom[amrlev][0].Domain());
124
125 for (MFIter mfi(a_psi_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
126 {
127 Box bx = mfi.growntilebox() & domain;
128
129 amrex::Array4<Set::Scalar> const& m_psi = (*(m_psi_mf[amrlev][0])).array(mfi);
130 amrex::Array4<const Set::Scalar> const& a_psi = a_psi_mf.array(mfi);
131
132 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
133 m_psi(i, j, k) = a_psi(i, j, k);
134 });
135 }
136 m_psi_set = true;
137}
138
139template<int SYM>
140void
141Elastic<SYM>::Fapply(int amrlev, int mglev, MultiFab& a_f, const MultiFab& a_u) const
142{
143 BL_PROFILE("Operator::Elastic::Fapply()");
144
145 amrex::Box domain(m_geom[amrlev][mglev].growPeriodicDomain(1));
146 domain.convert(amrex::IntVect::TheNodeVector());
147
148 amrex::Box stencilbox(m_geom[amrlev][mglev].growPeriodicDomain(2));
149 stencilbox.convert(amrex::IntVect::TheNodeVector());
150
151 const Real* DX = m_geom[amrlev][mglev].CellSize();
152
153 for (MFIter mfi(a_f, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
154 {
155 Box bx = mfi.validbox().grow(1) & domain;
156 amrex::Box tilebox = mfi.grownnodaltilebox() & bx;
157
158 amrex::Array4<MATRIX4> const& DDW = (*(m_ddw_mf[amrlev][mglev])).array(mfi);
159 amrex::Array4<const amrex::Real> const& U = a_u.array(mfi);
160 amrex::Array4<amrex::Real> const& F = a_f.array(mfi);
161 amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][mglev]->array(mfi);
162
163 const Dim3 lo = amrex::lbound(stencilbox), hi = amrex::ubound(stencilbox);
164
165 amrex::ParallelFor(tilebox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
166
167 Set::Vector f = Set::Vector::Zero();
168
169 Set::Vector u;
170 for (int p = 0; p < AMREX_SPACEDIM; p++) u(p) = U(i, j, k, p);
171
172
173 bool AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z)),
174 AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
175
176 // Determine if a special stencil will be necessary for first derivatives
177 std::array<Numeric::StencilType, AMREX_SPACEDIM>
178 sten = Numeric::GetStencil(i, j, k, stencilbox);
179
180 // The displacement gradient tensor
181 Set::Matrix gradu; // gradu(i,j) = u_{i,j)
182
183 // Fill gradu
184 for (int p = 0; p < AMREX_SPACEDIM; p++)
185 {
186 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(U, i, j, k, p, DX, sten));,
187 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(U, i, j, k, p, DX, sten));,
188 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(U, i, j, k, p, DX, sten)););
189 }
190
191 Set::Scalar psi_avg = 1.0;
192 if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
193
194 // Stress tensor computed using the model fab
195 Set::Matrix sig = (DDW(i, j, k) * gradu) * psi_avg;
196
197 // Boundary conditions
198 /// \todo Important: we need a way to handle corners and edges.
199 amrex::IntVect m(AMREX_D_DECL(i, j, k));
200 if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
201 {
202 f = (*m_bc)(u, gradu, sig, i, j, k, stencilbox);
203 }
204 else
205 {
206
207
208 // The gradient of the displacement gradient tensor
209 // TODO - replace with this call. But not for this PR
210 //Set::Matrix3 gradgradu = Numeric::Hessian(U,i,j,k,DX,sten); // gradgradu[k](l,j) = u_{k,lj}
211 Set::Matrix3 gradgradu; // gradgradu[k](l,j) = u_{k,lj}
212
213 // Fill gradu and gradgradu
214 for (int p = 0; p < AMREX_SPACEDIM; p++)
215 {
216 // Diagonal terms:
217 AMREX_D_TERM(
218 gradgradu(p, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(U, i, j, k, p, DX));,
219 gradgradu(p, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(U, i, j, k, p, DX));,
220 gradgradu(p, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(U, i, j, k, p, DX)););
221
222 // Off-diagonal terms:
223 AMREX_D_TERM(
224 ,// 2D
225 gradgradu(p, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(U, i, j, k, p, DX));
226 gradgradu(p, 1, 0) = gradgradu(p, 0, 1);
227 ,// 3D
228 gradgradu(p, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(U, i, j, k, p, DX));
229 gradgradu(p, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(U, i, j, k, p, DX));
230 gradgradu(p, 2, 0) = gradgradu(p, 0, 2);
231 gradgradu(p, 2, 1) = gradgradu(p, 1, 2););
232 }
233
234 //
235 // Operator
236 //
237 // The return value is
238 // f = C(grad grad u) + grad(C)*grad(u)
239 // In index notation
240 // f_i = C_{ijkl,j} u_{k,l} + C_{ijkl}u_{k,lj}
241 //
242
243 f = (DDW(i, j, k) * gradgradu) * psi_avg;
244
245 if (!m_uniform)
246 {
247 MATRIX4
248 AMREX_D_DECL(Cgrad1 = (Numeric::Stencil<MATRIX4, 1, 0, 0>::D(DDW, i, j, k, 0, DX, sten)),
249 Cgrad2 = (Numeric::Stencil<MATRIX4, 0, 1, 0>::D(DDW, i, j, k, 0, DX, sten)),
250 Cgrad3 = (Numeric::Stencil<MATRIX4, 0, 0, 1>::D(DDW, i, j, k, 0, DX, sten)));
251 f += (AMREX_D_TERM((Cgrad1 * gradu).col(0),
252 +(Cgrad2 * gradu).col(1),
253 +(Cgrad3 * gradu).col(2))) * (psi_avg);
254 }
255 if (m_psi_set)
256 {
257 Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, DX);
258 gradpsi *= (1.0 - m_psi_small);
259 f += (DDW(i, j, k) * gradu) * gradpsi;
260 }
261
262 if (std::isnan(f(0)) || std::isnan(f(1)))
263 {
264 Util::Message(INFO," ================= ");
265 Util::Message(INFO,"amrlev=",amrlev);
266 Util::Message(INFO,"mglev=",mglev);
267 Util::Message(INFO,"i=",i," j=",j);
268 Util::Message(INFO,"f: ",f.transpose());
269 Util::Message(INFO,"U(i,j,k): ",U(i,j,k,0)," ",U(i,j,k,1));
270 Util::Message(INFO,"U(i-1,j,k): ",U(i-1,j,k,0)," ",U(i-1,j,k,1));
271 Util::Message(INFO,"U(i+1,j,k): ",U(i+1,j,k,0)," ",U(i-1,j,k,1));
272 Util::Message(INFO,"U(i,j-1,k): ",U(i,j-1,k,0)," ",U(i,j-1,k,1));
273 Util::Message(INFO,"U(i,j+1,k): ",U(i,j+1,k,0)," ",U(i,j+1,k,1));
274 Util::Message(INFO,"gradu: ",gradu);
275 Util::Message(INFO,"gradgradu[0]: ",gradgradu[0]);
276 Util::Message(INFO,"gradgradu[1]: ",gradgradu[1]);
277 Util::Message(INFO,"DDW (i ,j ): ",DDW(i,j,k));
278 Util::Message(INFO,"DDW (i-1,j ): ",DDW(i-1,j,k));
279 Util::Message(INFO,"DDW (i+1,j ): ",DDW(i+1,j,k));
280 Util::Message(INFO,"DDW (i ,j+1): ",DDW(i,j+1,k));
281 Util::Message(INFO,"DDW (i ,j-1): ",DDW(i,j-1,k));
282 Util::Message(INFO,"psi_av: ",psi_avg);
283 Util::Message(INFO,"psi_set: ",m_psi_set);
284 Util::Message(INFO," ================= ");
285 Util::Abort(INFO);
286 }
287
288 }
289 AMREX_D_TERM(F(i, j, k, 0) = f[0];, F(i, j, k, 1) = f[1];, F(i, j, k, 2) = f[2];);
290 });
291 }
292}
293
294
295
296template<int SYM>
297void
298Elastic<SYM>::Diagonal(int amrlev, int mglev, MultiFab& a_diag)
299{
300 BL_PROFILE("Operator::Elastic::Diagonal()");
301
302 amrex::Box domain(m_geom[amrlev][mglev].growPeriodicDomain(1));
303 domain.convert(amrex::IntVect::TheNodeVector());
304
305 amrex::Box stencilbox(m_geom[amrlev][mglev].growPeriodicDomain(2));
306 stencilbox.convert(amrex::IntVect::TheNodeVector());
307
308 const Real* DX = m_geom[amrlev][mglev].CellSize();
309
310 for (MFIter mfi(a_diag, false); mfi.isValid(); ++mfi)
311 {
312 Box bx = mfi.validbox().grow(1) & domain;
313 amrex::Box tilebox = mfi.grownnodaltilebox() & bx;
314
315 amrex::Array4<MATRIX4> const& DDW = (*(m_ddw_mf[amrlev][mglev])).array(mfi);
316 amrex::Array4<Set::Scalar> const& diag = a_diag.array(mfi);
317 amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][mglev]->array(mfi);
318
319 const Dim3 lo = amrex::lbound(stencilbox), hi = amrex::ubound(stencilbox);
320
321 amrex::ParallelFor(tilebox, [=] AMREX_GPU_DEVICE(int i, int j, int k)
322 {
323
324 // Determine if a special stencil will be necessary for first derivatives
325 std::array<Numeric::StencilType, AMREX_SPACEDIM>
326 sten = Numeric::GetStencil(i, j, k, stencilbox);
327
328 // gradu(i,j) = u_{i,j)
329 std::array<Set::Matrix,AMREX_SPACEDIM> gradu = Numeric::Gradient_Diagonal<Set::Matrix>(DX, sten);
330
331 // gradgradu[k](l,j) = u_{k,lj}
332 std::array<Set::Matrix3,AMREX_SPACEDIM> gradgradu = Numeric::Gradient_Diagonal<Set::Matrix3>(DX);
333
334
335 Set::Vector f = Set::Vector::Zero();
336
337 bool
338 AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z)),
339 AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
340
341 Set::Scalar psi_avg = 1.0;
342 if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
343
344
345
346 for (int p = 0; p < AMREX_SPACEDIM; p++)
347 {
348
349 diag(i, j, k, p) = 0.0;
350
351
352 amrex::IntVect m(AMREX_D_DECL(i, j, k));
353 if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
354 {
355 Set::Matrix sig = DDW(i, j, k) * gradu[p] * psi_avg;
356 Set::Vector u = Set::Vector::Zero();
357 u(p) = 1.0;
358 f = (*m_bc)(u, gradu[p], sig, i, j, k, stencilbox);
359 diag(i, j, k, p) = f(p);
360 }
361 else
362 {
363 Set::Vector f = (DDW(i, j, k) * gradgradu[p]) * psi_avg;
364 diag(i, j, k, p) += f(p);
365 }
366
367#ifdef AMREX_DEBUG
368 if (std::isnan(diag(i, j, k, p))) Util::Abort(INFO, "diagonal is nan at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
369 if (std::isinf(diag(i, j, k, p))) Util::Abort(INFO, "diagonal is inf at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
370 if (diag(i, j, k, p) == 0) Util::Abort(INFO, "diagonal is zero at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
371#endif
372
373 }
374 });
375 }
376
377 a_diag.FillBoundaryAndSync(Geom(amrlev,mglev).periodicity());
378 nodalSync(amrlev,mglev,a_diag);
379}
380
381
382template<int SYM>
383void
384Elastic<SYM>::Error0x(int amrlev, int mglev, MultiFab& R0x, const MultiFab& x) const
385{
386 BL_PROFILE("Operator::Elastic::Error0x()");
387 Util::Message(INFO);
388
389 int ncomp = x.nComp();//getNComp();
390 int nghost = x.nGrow();
391
392 if (!m_diagonal_computed)
393 Util::Abort(INFO, "Operator::Diagonal() must be called before using normalize");
394
395 amrex::MultiFab D0x(x.boxArray(), x.DistributionMap(), ncomp, nghost);
396 amrex::MultiFab AD0x(x.boxArray(), x.DistributionMap(), ncomp, nghost);
397
398 amrex::MultiFab::Copy(D0x, x, 0, 0, ncomp, nghost); // D0x = x
399 amrex::MultiFab::Divide(D0x, *m_diag[amrlev][mglev], 0, 0, ncomp, 0); // D0x = x/diag
400 amrex::MultiFab::Copy(AD0x, D0x, 0, 0, ncomp, nghost); // AD0x = D0x
401
402 Fapply(amrlev, mglev, AD0x, D0x); // AD0x = A * D0 * x
403
404 amrex::MultiFab::Copy(R0x, x, 0, 0, ncomp, nghost); // R0x = x
405 amrex::MultiFab::Subtract(R0x, AD0x, 0, 0, ncomp, nghost); // R0x = x - AD0x
406}
407
408
409template<int SYM>
410void
411Elastic<SYM>::FFlux(int /*amrlev*/, const MFIter& /*mfi*/,
412 const std::array<FArrayBox*, AMREX_SPACEDIM>& sigmafab,
413 const FArrayBox& /*ufab*/, const int /*face_only*/) const
414{
415 BL_PROFILE("Operator::Elastic::FFlux()");
416 Util::Message(INFO);
417 amrex::BaseFab<amrex::Real> AMREX_D_DECL(&fxfab = *sigmafab[0],
418 &fyfab = *sigmafab[1],
419 &fzfab = *sigmafab[2]);
420 AMREX_D_TERM(fxfab.setVal(0.0);,
421 fyfab.setVal(0.0);,
422 fzfab.setVal(0.0););
423
424}
425
426template<int SYM>
427void
428Elastic<SYM>::Strain(int amrlev,
429 amrex::MultiFab& a_eps,
430 const amrex::MultiFab& a_u,
431 bool voigt) const
432{
433 BL_PROFILE("Operator::Elastic::Strain()");
434
435 const amrex::Real* DX = m_geom[amrlev][0].CellSize();
436 amrex::Box domain(m_geom[amrlev][0].Domain());
437 domain.convert(amrex::IntVect::TheNodeVector());
438
439
440 for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
441 {
442 const Box& bx = mfi.tilebox();
443 amrex::Array4<amrex::Real> const& epsilon = a_eps.array(mfi);
444 amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
445 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
446 {
447 Set::Matrix gradu;
448
449 std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
450 = Numeric::GetStencil(i, j, k, domain);
451
452 // Fill gradu
453 for (int p = 0; p < AMREX_SPACEDIM; p++)
454 {
455 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
456 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
457 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
458 }
459
460 Set::Matrix eps = 0.5 * (gradu + gradu.transpose());
461
462 if (voigt)
463 {
464 AMREX_D_PICK(epsilon(i, j, k, 0) = eps(0, 0);
465 ,
466 epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(1, 1); epsilon(i, j, k, 2) = eps(0, 1);
467 ,
468 epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(1, 1); epsilon(i, j, k, 2) = eps(2, 2);
469 epsilon(i, j, k, 3) = eps(1, 2); epsilon(i, j, k, 4) = eps(2, 0); epsilon(i, j, k, 5) = eps(0, 1););
470 }
471 else
472 {
473 AMREX_D_PICK(epsilon(i, j, k, 0) = eps(0, 0);
474 ,
475 epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(0, 1);
476 epsilon(i, j, k, 2) = eps(1, 0); epsilon(i, j, k, 3) = eps(1, 1);
477 ,
478 epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(0, 1); epsilon(i, j, k, 2) = eps(0, 2);
479 epsilon(i, j, k, 3) = eps(1, 0); epsilon(i, j, k, 4) = eps(1, 1); epsilon(i, j, k, 5) = eps(1, 2);
480 epsilon(i, j, k, 6) = eps(2, 0); epsilon(i, j, k, 7) = eps(2, 1); epsilon(i, j, k, 8) = eps(2, 2););
481 }
482 });
483 }
484}
485
486
487template<int SYM>
488void
489Elastic<SYM>::Stress(int amrlev,
490 amrex::MultiFab& a_sigma,
491 const amrex::MultiFab& a_u,
492 bool voigt, bool a_homogeneous)
493{
494 BL_PROFILE("Operator::Elastic::Stress()");
495 SetHomogeneous(a_homogeneous);
496
497 const amrex::Real* DX = m_geom[amrlev][0].CellSize();
498 amrex::Box domain(m_geom[amrlev][0].Domain());
499 domain.convert(amrex::IntVect::TheNodeVector());
500
501 for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
502 {
503 const Box& bx = mfi.tilebox();
504 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& DDW = (*(m_ddw_mf[amrlev][0])).array(mfi);
505 amrex::Array4<amrex::Real> const& sigma = a_sigma.array(mfi);
506 amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][0]->array(mfi);
507 amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
508 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
509 {
510 Set::Matrix gradu;
511
512 std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
513 = Numeric::GetStencil(i, j, k, domain);
514
515 // Fill gradu
516 for (int p = 0; p < AMREX_SPACEDIM; p++)
517 {
518 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
519 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
520 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
521 }
522
523 Set::Scalar psi_avg = 1.0;
524 if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
525 Set::Matrix sig = (DDW(i, j, k) * gradu) * psi_avg;
526
527 if (voigt)
528 {
529 AMREX_D_PICK(sigma(i, j, k, 0) = sig(0, 0);
530 ,
531 sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(1, 1); sigma(i, j, k, 2) = sig(0, 1);
532 ,
533 sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(1, 1); sigma(i, j, k, 2) = sig(2, 2);
534 sigma(i, j, k, 3) = sig(1, 2); sigma(i, j, k, 4) = sig(2, 0); sigma(i, j, k, 5) = sig(0, 1););
535 }
536 else
537 {
538 AMREX_D_PICK(sigma(i, j, k, 0) = sig(0, 0);
539 ,
540 sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(0, 1);
541 sigma(i, j, k, 2) = sig(1, 0); sigma(i, j, k, 3) = sig(1, 1);
542 ,
543 sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(0, 1); sigma(i, j, k, 2) = sig(0, 2);
544 sigma(i, j, k, 3) = sig(1, 0); sigma(i, j, k, 4) = sig(1, 1); sigma(i, j, k, 5) = sig(1, 2);
545 sigma(i, j, k, 6) = sig(2, 0); sigma(i, j, k, 7) = sig(2, 1); sigma(i, j, k, 8) = sig(2, 2););
546 }
547 });
548 }
549}
550
551
552template<int SYM>
553void
554Elastic<SYM>::Energy(int amrlev,
555 amrex::MultiFab& a_energy,
556 const amrex::MultiFab& a_u, bool a_homogeneous)
557{
558 BL_PROFILE("Operator::Elastic::Energy()");
559 SetHomogeneous(a_homogeneous);
560
561 amrex::Box domain(m_geom[amrlev][0].Domain());
562 domain.convert(amrex::IntVect::TheNodeVector());
563
564 const amrex::Real* DX = m_geom[amrlev][0].CellSize();
565
566 for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
567 {
568 const Box& bx = mfi.tilebox();
569 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& DDW = (*(m_ddw_mf[amrlev][0])).array(mfi);
570 amrex::Array4<amrex::Real> const& energy = a_energy.array(mfi);
571 amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
572 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
573 {
574 Set::Matrix gradu;
575
576 std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
577 = Numeric::GetStencil(i, j, k, domain);
578
579 // Fill gradu
580 for (int p = 0; p < AMREX_SPACEDIM; p++)
581 {
582 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
583 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
584 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
585 }
586
587 Set::Matrix eps = .5 * (gradu + gradu.transpose());
588 Set::Matrix sig = DDW(i, j, k) * gradu;
589
590 // energy(i,j,k) = (gradu.transpose() * sig).trace();
591
592 //Util::Abort(INFO,"Fix this"); //
593 //energy(i,j,k) = C(i,j,k).W(gradu);
594 for (int m = 0; m < AMREX_SPACEDIM; m++)
595 {
596 for (int n = 0; n < AMREX_SPACEDIM; n++)
597 {
598 energy(i, j, k) += .5 * sig(m, n) * eps(m, n);
599 }
600 }
601 });
602 }
603}
604
605template<int SYM>
606void
607Elastic<SYM>::averageDownCoeffs()
608{
609 BL_PROFILE("Elastic::averageDownCoeffs()");
610
611 if (m_average_down_coeffs)
612 for (int amrlev = m_num_amr_levels - 1; amrlev > 0; --amrlev)
613 averageDownCoeffsDifferentAmrLevels(amrlev);
614
615 averageDownCoeffsSameAmrLevel(0);
616 for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
617 {
618 for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
619 {
620 if (m_ddw_mf[amrlev][mglev]) {
621 FillBoundaryCoeff(*m_ddw_mf[amrlev][mglev], Geom(amrlev,mglev).periodicity());
622 FillBoundaryCoeff(*m_psi_mf[amrlev][mglev], Geom(amrlev,mglev).periodicity());
623 }
624 }
625 }
626}
627
628template<int SYM>
629void
630Elastic<SYM>::averageDownCoeffsDifferentAmrLevels(int fine_amrlev)
631{
632 BL_PROFILE("Operator::Elastic::averageDownCoeffsDifferentAmrLevels()");
633 Util::Assert(INFO, TEST(fine_amrlev > 0));
634
635 const int crse_amrlev = fine_amrlev - 1;
636 const int ncomp = 1;
637
638 MultiTab& crse_ddw = *m_ddw_mf[crse_amrlev][0];
639 MultiTab& fine_ddw = *m_ddw_mf[fine_amrlev][0];
640
641 amrex::Box cdomain(m_geom[crse_amrlev][0].Domain());
642 cdomain.convert(amrex::IntVect::TheNodeVector());
643
644 const Geometry& cgeom = m_geom[crse_amrlev][0];
645
646 const BoxArray& fba = fine_ddw.boxArray();
647 const DistributionMapping& fdm = fine_ddw.DistributionMap();
648
649 MultiTab fine_ddw_for_coarse(amrex::coarsen(fba, 2), fdm, ncomp, 2);
650 fine_ddw_for_coarse.ParallelCopy(crse_ddw, 0, 0, ncomp, 0, 0, cgeom.periodicity());
651
652 const int coarse_fine_node = 1;
653 const int fine_fine_node = 2;
654
655 amrex::iMultiFab nodemask(amrex::coarsen(fba, 2), fdm, 1, 2);
656 nodemask.ParallelCopy(*m_nd_fine_mask[crse_amrlev], 0, 0, 1, 0, 0, cgeom.periodicity());
657
658 amrex::iMultiFab cellmask(amrex::convert(amrex::coarsen(fba, 2), amrex::IntVect::TheCellVector()), fdm, 1, 2);
659 cellmask.ParallelCopy(*m_cc_fine_mask[crse_amrlev], 0, 0, 1, 1, 1, cgeom.periodicity());
660
661 for (MFIter mfi(fine_ddw_for_coarse, false); mfi.isValid(); ++mfi)
662 {
663 const Box& bx = mfi.validbox();
664
665 amrex::Array4<const int> const& nmask = nodemask.array(mfi);
666 //amrex::Array4<const int> const& cmask = cellmask.array(mfi);
667
668 amrex::Array4<MATRIX4> const& cdata = fine_ddw_for_coarse.array(mfi);
669 amrex::Array4<const MATRIX4> const& fdata = fine_ddw.array(mfi);
670
671 const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
672
673 for (int n = 0; n < fine_ddw.nComp(); n++)
674 {
675 // I,J,K == coarse coordinates
676 // i,j,k == fine coordinates
677 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
678 int i = I * 2, j = J * 2, k = K * 2;
679
680 if (nmask(I, J, K) == fine_fine_node || nmask(I, J, K) == coarse_fine_node)
681 {
682 if ((I == lo.x || I == hi.x) &&
683 (J == lo.y || J == hi.y) &&
684 (K == lo.z || K == hi.z)) // Corner
685 cdata(I, J, K, n) = fdata(i, j, k, n);
686 else if ((J == lo.y || J == hi.y) &&
687 (K == lo.z || K == hi.z)) // X edge
688 cdata(I, J, K, n) = fdata(i - 1, j, k, n) * 0.25 + fdata(i, j, k, n) * 0.5 + fdata(i + 1, j, k, n) * 0.25;
689 else if ((K == lo.z || K == hi.z) &&
690 (I == lo.x || I == hi.x)) // Y edge
691 cdata(I, J, K, n) = fdata(i, j - 1, k, n) * 0.25 + fdata(i, j, k, n) * 0.5 + fdata(i, j + 1, k, n) * 0.25;
692 else if ((I == lo.x || I == hi.x) &&
693 (J == lo.y || J == hi.y)) // Z edge
694 cdata(I, J, K, n) = fdata(i, j, k - 1, n) * 0.25 + fdata(i, j, k, n) * 0.5 + fdata(i, j, k + 1, n) * 0.25;
695 else if (I == lo.x || I == hi.x) // X face
696 cdata(I, J, K, n) =
697 (fdata(i, j - 1, k - 1, n) + fdata(i, j, k - 1, n) * 2.0 + fdata(i, j + 1, k - 1, n)
698 + fdata(i, j - 1, k, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i, j + 1, k, n) * 2.0
699 + fdata(i, j - 1, k + 1, n) + fdata(i, j, k + 1, n) * 2.0 + fdata(i, j + 1, k + 1, n)) / 16.0;
700 else if (J == lo.y || J == hi.y) // Y face
701 cdata(I, J, K, n) =
702 (fdata(i - 1, j, k - 1, n) + fdata(i - 1, j, k, n) * 2.0 + fdata(i - 1, j, k + 1, n)
703 + fdata(i, j, k - 1, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i, j, k + 1, n) * 2.0
704 + fdata(i + 1, j, k - 1, n) + fdata(i + 1, j, k, n) * 2.0 + fdata(i + 1, j, k + 1, n)) / 16.0;
705 else if (K == lo.z || K == hi.z) // Z face
706 cdata(I, J, K, n) =
707 (fdata(i - 1, j - 1, k, n) + fdata(i, j - 1, k, n) * 2.0 + fdata(i + 1, j - 1, k, n)
708 + fdata(i - 1, j, k, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i + 1, j, k, n) * 2.0
709 + fdata(i - 1, j + 1, k, n) + fdata(i, j + 1, k, n) * 2.0 + fdata(i + 1, j + 1, k, n)) / 16.0;
710 else // Interior
711 cdata(I, J, K, n) =
712 (fdata(i - 1, j - 1, k - 1, n) + fdata(i - 1, j - 1, k + 1, n) + fdata(i - 1, j + 1, k - 1, n) + fdata(i - 1, j + 1, k + 1, n) +
713 fdata(i + 1, j - 1, k - 1, n) + fdata(i + 1, j - 1, k + 1, n) + fdata(i + 1, j + 1, k - 1, n) + fdata(i + 1, j + 1, k + 1, n)) / 64.0
714 +
715 (fdata(i, j - 1, k - 1, n) + fdata(i, j - 1, k + 1, n) + fdata(i, j + 1, k - 1, n) + fdata(i, j + 1, k + 1, n) +
716 fdata(i - 1, j, k - 1, n) + fdata(i + 1, j, k - 1, n) + fdata(i - 1, j, k + 1, n) + fdata(i + 1, j, k + 1, n) +
717 fdata(i - 1, j - 1, k, n) + fdata(i - 1, j + 1, k, n) + fdata(i + 1, j - 1, k, n) + fdata(i + 1, j + 1, k, n)) / 32.0
718 +
719 (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
720 fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
721 +
722 fdata(i, j, k, n) / 8.0;
723
724#ifdef AMREX_DEBUG
725 if (cdata(I, J, K).contains_nan()) Util::Abort(INFO, "restricted model is nan at (", i, ",", j, ",", k, "), fine_amrlev=", fine_amrlev);
726#endif
727 }
728
729 });
730 }
731 }
732
733 // Copy the fine residual restricted onto the coarse grid
734 // into the final residual.
735
736 crse_ddw.ParallelCopy(fine_ddw_for_coarse, 0, 0, ncomp, 0, 0, cgeom.periodicity());
737 //const int mglev = 0;
738 //Util::RealFillBoundary(crse_ddw, m_geom[crse_amrlev][mglev]);
739
740 FillBoundaryCoeff(crse_ddw,Geom(fine_amrlev,0).periodicity());
741
742}
743
744
745
746template<int SYM>
747void
748Elastic<SYM>::averageDownCoeffsSameAmrLevel(int amrlev)
749{
750 BL_PROFILE("Elastic::averageDownCoeffsSameAmrLevel()");
751
752 for (int mglev = 1; mglev < m_num_mg_levels[amrlev]; ++mglev)
753 {
754 amrex::Box cdomain(m_geom[amrlev][mglev].growPeriodicDomain(2));
755 cdomain.convert(amrex::IntVect::TheNodeVector());
756 amrex::Box fdomain(m_geom[amrlev][mglev - 1].Domain());
757 fdomain.convert(amrex::IntVect::TheNodeVector());
758
759 MultiTab& crse = *m_ddw_mf[amrlev][mglev];
760 MultiTab& fine = *m_ddw_mf[amrlev][mglev - 1];
761
762 amrex::BoxArray crseba = crse.boxArray();
763 amrex::BoxArray fineba = fine.boxArray();
764
765 BoxArray newba = crseba;
766 newba.refine(2);
767 MultiTab fine_on_crseba;
768 fine_on_crseba.define(newba, crse.DistributionMap(), 1, 4);
769 fine_on_crseba.ParallelCopy(fine, 0, 0, 1, 2, 4, m_geom[amrlev][mglev-1].periodicity());
770 /* ine_on_crseba.FillBoundaryAndSync(m_geom[amrlev][mglev-1].periodicity()); */
771
772 for (MFIter mfi(crse, false); mfi.isValid(); ++mfi)
773 {
774
775 Box bx = mfi.grownnodaltilebox() & cdomain;
776 /*Box bx = mfi.grownnodaltilebox(-1,1) & cdomain;*/
777
778 amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, SYM>> const& fdata = fine_on_crseba.array(mfi);
779 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& cdata = crse.array(mfi);
780
781 const Dim3 lo = amrex::lbound(bx), hi = amrex::ubound(bx);
782 /*const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);*/
783
784 // I,J,K == coarse coordinates
785 // i,j,k == fine coordinates
786 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
787 int i = 2 * I, j = 2 * J, k = 2 * K;
788
789 if ((I == lo.x || I == hi.x) &&
790 (J == lo.y || J == hi.y) &&
791 (K == lo.z || K == hi.z)) // Corner
792 cdata(I, J, K) = fdata(i, j, k);
793 else if ((J == lo.y || J == hi.y) &&
794 (K == lo.z || K == hi.z)) // X edge
795 cdata(I, J, K) = fdata(i - 1, j, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i + 1, j, k) * 0.25;
796 else if ((K == lo.z || K == hi.z) &&
797 (I == lo.x || I == hi.x)) // Y edge
798 cdata(I, J, K) = fdata(i, j - 1, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j + 1, k) * 0.25;
799 else if ((I == lo.x || I == hi.x) &&
800 (J == lo.y || J == hi.y)) // Z edge
801 cdata(I, J, K) = fdata(i, j, k - 1) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j, k + 1) * 0.25;
802 else if (I == lo.x || I == hi.x) // X face
803 cdata(I, J, K) =
804 (fdata(i, j - 1, k - 1) + fdata(i, j, k - 1) * 2.0 + fdata(i, j + 1, k - 1)
805 + fdata(i, j - 1, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j + 1, k) * 2.0
806 + fdata(i, j - 1, k + 1) + fdata(i, j, k + 1) * 2.0 + fdata(i, j + 1, k + 1)) / 16.0;
807 else if (J == lo.y || J == hi.y) // Y face
808 cdata(I, J, K) =
809 (fdata(i - 1, j, k - 1) + fdata(i - 1, j, k) * 2.0 + fdata(i - 1, j, k + 1)
810 + fdata(i, j, k - 1) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j, k + 1) * 2.0
811 + fdata(i + 1, j, k - 1) + fdata(i + 1, j, k) * 2.0 + fdata(i + 1, j, k + 1)) / 16.0;
812 else if (K == lo.z || K == hi.z) // Z face
813 cdata(I, J, K) =
814 (fdata(i - 1, j - 1, k) + fdata(i, j - 1, k) * 2.0 + fdata(i + 1, j - 1, k)
815 + fdata(i - 1, j, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i + 1, j, k) * 2.0
816 + fdata(i - 1, j + 1, k) + fdata(i, j + 1, k) * 2.0 + fdata(i + 1, j + 1, k)) / 16.0;
817 else // Interior
818 cdata(I, J, K) =
819 (fdata(i - 1, j - 1, k - 1) + fdata(i - 1, j - 1, k + 1) + fdata(i - 1, j + 1, k - 1) + fdata(i - 1, j + 1, k + 1) +
820 fdata(i + 1, j - 1, k - 1) + fdata(i + 1, j - 1, k + 1) + fdata(i + 1, j + 1, k - 1) + fdata(i + 1, j + 1, k + 1)) / 64.0
821 +
822 (fdata(i, j - 1, k - 1) + fdata(i, j - 1, k + 1) + fdata(i, j + 1, k - 1) + fdata(i, j + 1, k + 1) +
823 fdata(i - 1, j, k - 1) + fdata(i + 1, j, k - 1) + fdata(i - 1, j, k + 1) + fdata(i + 1, j, k + 1) +
824 fdata(i - 1, j - 1, k) + fdata(i - 1, j + 1, k) + fdata(i + 1, j - 1, k) + fdata(i + 1, j + 1, k)) / 32.0
825 +
826 (fdata(i - 1, j, k) + fdata(i, j - 1, k) + fdata(i, j, k - 1) +
827 fdata(i + 1, j, k) + fdata(i, j + 1, k) + fdata(i, j, k + 1)) / 16.0
828 +
829 fdata(i, j, k) / 8.0;
830
831#ifdef AMREX_DEBUG
832 if (cdata(I, J, K).contains_nan()) Util::Abort(INFO, "restricted model is nan at crse coordinates (I=", I, ",J=", J, ",K=", k, "), amrlev=", amrlev, " interpolating from mglev", mglev - 1, " to ", mglev);
833#endif
834 });
835 }
836 FillBoundaryCoeff(crse, Geom(amrlev,mglev).periodicity());
837
838
839 if (!m_psi_set) continue;
840
841 amrex::Box cdomain_cell(m_geom[amrlev][mglev].Domain());
842 amrex::Box fdomain_cell(m_geom[amrlev][mglev - 1].Domain());
843 MultiFab& crse_psi = *m_psi_mf[amrlev][mglev];
844 MultiFab& fine_psi = *m_psi_mf[amrlev][mglev - 1];
845 MultiFab fine_psi_on_crseba;
846 fine_psi_on_crseba.define(newba.convert(amrex::IntVect::TheCellVector()), crse_psi.DistributionMap(), 1, 1);
847 fine_psi_on_crseba.ParallelCopy(fine_psi, 0, 0, 1, 1, 1, m_geom[amrlev][mglev].periodicity());
848
849 for (MFIter mfi(crse_psi, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
850 {
851 Box bx = mfi.tilebox();
852 bx = bx & cdomain_cell;
853
854 amrex::Array4<const Set::Scalar> const& fdata = fine_psi_on_crseba.array(mfi);
855 amrex::Array4<Set::Scalar> const& cdata = crse_psi.array(mfi);
856
857 const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
858
859 // I,J,K == coarse coordinates
860 // i,j,k == fine coordinates
861 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
862 int i = 2 * I, j = 2 * J, k = 2 * K;
863
864 if ((I == lo.x || I == hi.x) &&
865 (J == lo.y || J == hi.y) &&
866 (K == lo.z || K == hi.z)) // Corner
867 cdata(I, J, K) = fdata(i, j, k);
868 else if ((J == lo.y || J == hi.y) &&
869 (K == lo.z || K == hi.z)) // X edge
870 cdata(I, J, K) = fdata(i - 1, j, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i + 1, j, k) * 0.25;
871 else if ((K == lo.z || K == hi.z) &&
872 (I == lo.x || I == hi.x)) // Y edge
873 cdata(I, J, K) = fdata(i, j - 1, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j + 1, k) * 0.25;
874 else if ((I == lo.x || I == hi.x) &&
875 (J == lo.y || J == hi.y)) // Z edge
876 cdata(I, J, K) = fdata(i, j, k - 1) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j, k + 1) * 0.25;
877 else if (I == lo.x || I == hi.x) // X face
878 cdata(I, J, K) =
879 (fdata(i, j - 1, k - 1) + fdata(i, j, k - 1) * 2.0 + fdata(i, j + 1, k - 1)
880 + fdata(i, j - 1, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j + 1, k) * 2.0
881 + fdata(i, j - 1, k + 1) + fdata(i, j, k + 1) * 2.0 + fdata(i, j + 1, k + 1)) / 16.0;
882 else if (J == lo.y || J == hi.y) // Y face
883 cdata(I, J, K) =
884 (fdata(i - 1, j, k - 1) + fdata(i - 1, j, k) * 2.0 + fdata(i - 1, j, k + 1)
885 + fdata(i, j, k - 1) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j, k + 1) * 2.0
886 + fdata(i + 1, j, k - 1) + fdata(i + 1, j, k) * 2.0 + fdata(i + 1, j, k + 1)) / 16.0;
887 else if (K == lo.z || K == hi.z) // Z face
888 cdata(I, J, K) =
889 (fdata(i - 1, j - 1, k) + fdata(i, j - 1, k) * 2.0 + fdata(i + 1, j - 1, k)
890 + fdata(i - 1, j, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i + 1, j, k) * 2.0
891 + fdata(i - 1, j + 1, k) + fdata(i, j + 1, k) * 2.0 + fdata(i + 1, j + 1, k)) / 16.0;
892 else // Interior
893 cdata(I, J, K) =
894 (fdata(i - 1, j - 1, k - 1) + fdata(i - 1, j - 1, k + 1) + fdata(i - 1, j + 1, k - 1) + fdata(i - 1, j + 1, k + 1) +
895 fdata(i + 1, j - 1, k - 1) + fdata(i + 1, j - 1, k + 1) + fdata(i + 1, j + 1, k - 1) + fdata(i + 1, j + 1, k + 1)) / 64.0
896 +
897 (fdata(i, j - 1, k - 1) + fdata(i, j - 1, k + 1) + fdata(i, j + 1, k - 1) + fdata(i, j + 1, k + 1) +
898 fdata(i - 1, j, k - 1) + fdata(i + 1, j, k - 1) + fdata(i - 1, j, k + 1) + fdata(i + 1, j, k + 1) +
899 fdata(i - 1, j - 1, k) + fdata(i - 1, j + 1, k) + fdata(i + 1, j - 1, k) + fdata(i + 1, j + 1, k)) / 32.0
900 +
901 (fdata(i - 1, j, k) + fdata(i, j - 1, k) + fdata(i, j, k - 1) +
902 fdata(i + 1, j, k) + fdata(i, j + 1, k) + fdata(i, j, k + 1)) / 16.0
903 +
904 fdata(i, j, k) / 8.0;
905 });
906 }
907 FillBoundaryCoeff(crse_psi, Geom(amrlev,mglev).periodicity());
908
909 }
910}
911
912template<int SYM>
913void
914Elastic<SYM>::FillBoundaryCoeff(MultiTab& sigma, const amrex::Periodicity& p)
915{
916 sigma.setMultiGhost(true);
917 sigma.FillBoundaryAndSync(p);
918}
919
920template<int SYM>
921void
922Elastic<SYM>::FillBoundaryCoeff(MultiFab& psi, const amrex::Periodicity& p)
923{
924 psi.setMultiGhost(true);
925 psi.FillBoundaryAndSync(p);
926}
927
928template class Elastic<Set::Sym::Major>;
929template class Elastic<Set::Sym::Isotropic>;
930template class Elastic<Set::Sym::MajorMinor>;
931template class Elastic<Set::Sym::Diagonal>;
932
933}
TEST
#define TEST(x)
Definition Util.H:25
INFO
#define INFO
Definition Util.H:24
Operator::Elastic::define
void define(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const LPInfo &a_info=LPInfo(), const Vector< FabFactory< FArrayBox > const * > &a_factory={})
Definition Elastic.cpp:26
Operator::Elastic::SetModel
void SetModel(Set::Matrix4< AMREX_SPACEDIM, SYM > &a_model)
Definition Elastic.cpp:59
Operator::Elastic< T::sym >::MultiTab
amrex::FabArray< TArrayBox > MultiTab
Definition Elastic.H:26
Operator::Elastic::Stress
void Stress(int amrlev, amrex::MultiFab &sigmafab, const amrex::MultiFab &ufab, bool voigt=false, bool a_homogeneous=false)
Compute stress given the displacement field by.
Definition Elastic.cpp:489
Operator::Elastic::FFlux
virtual void FFlux(int amrlev, const MFIter &mfi, const std::array< FArrayBox *, AMREX_SPACEDIM > &flux, const FArrayBox &sol, const int face_only=0) const final
Definition Elastic.cpp:411
Operator::Elastic::Error0x
void Error0x(int amrlev, int mglev, MultiFab &R0x, const MultiFab &x) const
Definition Elastic.cpp:384
Operator::Elastic::FillBoundaryCoeff
void FillBoundaryCoeff(MultiTab &sigma, const amrex::Periodicity &p)
Definition Elastic.cpp:914
Operator::Elastic::averageDownCoeffsSameAmrLevel
void averageDownCoeffsSameAmrLevel(int amrlev)
Update coarse-level AMR coefficients with data from fine level.
Definition Elastic.cpp:748
Operator::Elastic::~Elastic
virtual ~Elastic()
Definition Elastic.cpp:21
Operator::Elastic::Strain
void Strain(int amrlev, amrex::MultiFab &epsfab, const amrex::MultiFab &ufab, bool voigt=false) const
Compute strain given the displacement field by.
Definition Elastic.cpp:428
Operator::Elastic::Energy
void Energy(int amrlev, amrex::MultiFab &energy, const amrex::MultiFab &u, bool a_homogeneous=false)
Compute energy density given the displacement field by.
Definition Elastic.cpp:554
Operator::Elastic::averageDownCoeffsDifferentAmrLevels
void averageDownCoeffsDifferentAmrLevels(int fine_amrlev)
Definition Elastic.cpp:630
Operator::Elastic::Fapply
virtual void Fapply(int amrlev, int mglev, MultiFab &out, const MultiFab &in) const override final
Definition Elastic.cpp:141
Operator::Elastic::Diagonal
virtual void Diagonal(int amrlev, int mglev, amrex::MultiFab &diag) override
Definition Elastic.cpp:298
Operator::Elastic::averageDownCoeffs
virtual void averageDownCoeffs() override
Definition Elastic.cpp:607
Operator::Elastic::SetPsi
void SetPsi(int amrlev, const amrex::MultiFab &a_psi)
Definition Elastic.cpp:120
Numeric::GetStencil
static AMREX_FORCE_INLINE std::array< StencilType, AMREX_SPACEDIM > GetStencil(const int i, const int j, const int k, const amrex::Box domain)
Definition Stencil.H:45
Numeric::CellGradientOnNode
AMREX_FORCE_INLINE Set::Vector CellGradientOnNode(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Stencil.H:699
Operator
Documentation for operator namespace.
Definition Diagonal.cpp:14
Operator::AMREX_D_DECL
constexpr amrex::IntVect AMREX_D_DECL(Operator< Grid::Cell >::dx, Operator< Grid::Cell >::dy, Operator< Grid::Cell >::dz)
Set::Scalar
amrex::Real Scalar
Definition Base.H:18
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:19
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition Base.H:22
Util::Abort
void Abort(const char *msg)
Definition Util.cpp:268
Util::Assert
AMREX_FORCE_INLINE void Assert(std::string file, std::string func, int line, std::string smt, bool pass, Args const &... args)
Definition Util.H:58
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:129
Numeric::Interpolate::CellToNodeAverage
static AMREX_FORCE_INLINE T CellToNodeAverage(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition Stencil.H:1390
Generated by doxygen 1.9.8