Line data Source code
1 : //
2 : // Create a single laminate with specified orientation, thickness, and offset.
3 : //
4 :
5 : #ifndef IC_LAMINATE_H_
6 : #define IC_LAMINATE_H_
7 :
8 : #include "IC/IC.H"
9 : #include "Util/Util.H"
10 : #include "IO/ParmParse.H"
11 :
12 : /// \class Laminate
13 : /// \brief Initialize Laminates in a matrix
14 : namespace IC
15 : {
16 : class Laminate: public IC
17 : {
18 : public:
19 : static constexpr const char* name = "laminate";
20 : enum Mollifier { Dirac, Gaussian };
21 :
22 0 : Laminate(amrex::Vector<amrex::Geometry>& _geom): IC(_geom) {}
23 3 : Laminate(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp, std::string name): IC(_geom)
24 : {
25 3 : pp_queryclass(name, *this);
26 3 : }
27 :
28 39 : void Add(const int& lev, Set::Field<Set::Scalar>& a_field, Set::Scalar)
29 : {
30 105 : for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
31 : {
32 66 : amrex::Box bx = mfi.tilebox();
33 66 : bx.grow(a_field[lev]->nGrow());
34 66 : amrex::IndexType type = a_field[lev]->ixType();
35 :
36 66 : amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
37 66 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
38 :
39 73042 : Set::Vector x = Set::Position(i, j, k, geom[lev], type);
40 :
41 : // The base matrix is always field(i,j,k,0)
42 : // Inclusions start with field(i,j,k,1)
43 :
44 73042 : if (!singlefab)
45 : {
46 0 : Set::Scalar value = 0.0;
47 0 : for (int m = 0; m < number_of_inclusions; m++)
48 : {
49 0 : Set::Scalar t = std::abs((x - center[m]).transpose() * orientation[m]) - (thickness[m] / 2.0);
50 0 : field(i, j, k, m + 1) = 0.5 - 0.5 * std::erf(t / eps[m]);
51 0 : if (field(i, j, k, m + 1) < 0.) field(i, j, k, m + 1) = 0.;
52 0 : if (field(i, j, k, m + 1) > 1.) field(i, j, k, m + 1) = 1.;
53 0 : value += field(i, j, k, m + 1);
54 : }
55 0 : field(i, j, k, 0) = 1.0 - value;
56 0 : if (field(i, j, k, 0) < 0.) field(i, j, k, 0) = 0.;
57 0 : if (field(i, j, k, 0) > 1.) field(i, j, k, 0) = 1.;
58 : }
59 : else
60 : {
61 73042 : Set::Scalar t = std::abs((x - center[0]).transpose() * orientation[0]) - (thickness[0] / 2.0);
62 73042 : field(i, j, k, 0) = 0.5 - 0.5 * std::erf(t / eps[0]);
63 219126 : if (invert) field(i, j, k, 0) = 1.0 - field(i, j, k, 0);
64 : }
65 73042 : });
66 : }
67 39 : a_field[lev]->FillBoundary();
68 39 : };
69 :
70 : private:
71 : int number_of_inclusions = -1;
72 : amrex::Vector<Set::Vector> center;
73 : amrex::Vector<Set::Vector> orientation;
74 : amrex::Vector<Set::Vector> normal;
75 : amrex::Vector<Set::Scalar> eps;
76 : amrex::Vector<Set::Scalar> thickness;
77 : Mollifier moll = Mollifier::Gaussian;
78 : bool singlefab = false;
79 : bool invert = false;
80 :
81 : public:
82 3 : static void Parse(Laminate& value, IO::ParmParse& pp)
83 : {
84 : // How many laminates (MUST be greater than or equal to 1).
85 3 : pp_query_default("number_of_inclusions", value.number_of_inclusions, 1);
86 :
87 3 : if (value.number_of_inclusions < 1)
88 0 : Util::Abort(INFO, "Number of inclusions must be at least 1. Aborting.");
89 :
90 6 : amrex::Vector<Set::Scalar> a_center;
91 : // (x,y,[z]) values for the center point of the laminate
92 3 : if (pp.contains("center")) pp_queryarr("center", a_center);
93 :
94 3 : if (a_center.size() != value.number_of_inclusions * AMREX_SPACEDIM) value.center.push_back(Set::Vector::Zero());
95 : else
96 : {
97 0 : for (int i = 0; i < a_center.size(); i += AMREX_SPACEDIM)
98 0 : value.center.push_back(Set::Vector(AMREX_D_DECL(a_center[i], a_center[i + 1], a_center[i + 2])));
99 : }
100 :
101 6 : amrex::Vector<Set::Scalar> a_thickness;
102 : // thickness of the laminate
103 3 : if (pp.contains("thickness")) pp_queryarr("thickness", a_thickness);
104 :
105 3 : if (a_thickness.size() != value.number_of_inclusions && a_thickness.size() != 1)
106 0 : Util::Abort(INFO, "Thickness of each inclusion must be specified");
107 :
108 3 : if (a_thickness.size() == 1)
109 : {
110 3 : if (a_thickness[0] <= 0.0) Util::Abort(INFO, "Invalid value of inclusion thickness");
111 6 : for (int i = 0; i < value.number_of_inclusions; i++) value.thickness.push_back(a_thickness[0]);
112 : }
113 :
114 : else
115 : {
116 0 : for (int i = 0; i < value.number_of_inclusions; i++)
117 : {
118 0 : if (a_thickness[i] <= 0.0) Util::Abort(INFO, "Invalid value of inclusion ", i + 1, " thickness");
119 0 : value.thickness.push_back(a_thickness[i]);
120 : }
121 : }
122 :
123 6 : amrex::Vector<Set::Scalar> a_orientation;
124 : // Vector normal to the interface of the laminate
125 3 : pp_queryarr("orientation", a_orientation);
126 :
127 3 : if (a_orientation.size() != value.number_of_inclusions * AMREX_SPACEDIM && a_orientation.size() != AMREX_SPACEDIM)
128 0 : Util::Abort(INFO, "Orientation of each inclusion must be specified");
129 :
130 3 : if (a_orientation.size() == AMREX_SPACEDIM)
131 6 : for (int i = 0; i < value.number_of_inclusions; i++)
132 3 : value.orientation.push_back(Set::Vector(AMREX_D_DECL(a_orientation[0], a_orientation[1], a_orientation[2])));
133 :
134 : else
135 0 : for (int i = 0; i < a_orientation.size(); i += AMREX_SPACEDIM)
136 0 : value.orientation.push_back(Set::Vector(AMREX_D_DECL(a_orientation[i], a_orientation[i + 1], a_orientation[i + 2])));
137 :
138 6 : for (int i = 0; i < value.orientation.size(); i++)
139 3 : if (value.orientation[i].lpNorm<2>() <= 0.) value.orientation[i] = Set::Vector::Random();
140 :
141 6 : for (int i = 0; i < value.orientation.size(); i++)
142 3 : value.orientation[i] = value.orientation[i] / value.orientation[i].lpNorm<2>();
143 :
144 3 : value.normal.resize(value.number_of_inclusions);
145 6 : for (int i = 0; i < value.orientation.size(); i++)
146 : {
147 3 : value.normal[i] = Set::Vector::Zero();
148 3 : if (value.orientation[i](0) != 0.)
149 : {
150 0 : AMREX_D_TERM(value.normal[i](0) = 1.;, value.normal[i](1) = 1.;, value.normal[i](2) = 1.;);
151 0 : value.normal[i](0) = -(AMREX_D_TERM(0., +value.orientation[i](1), +value.orientation[i](2))) / value.orientation[i](0);
152 0 : value.normal[i] = value.normal[i] / value.normal[i].lpNorm<2>();
153 : }
154 3 : else if (value.orientation[i](1) != 0.)
155 : {
156 3 : AMREX_D_TERM(value.normal[i](0) = 1.;, value.normal[i](1) = 1.;, value.normal[i](2) = 1.;);
157 3 : value.normal[i](1) = -(AMREX_D_TERM(value.orientation[i](0), +0.0, +value.orientation[i](2))) / value.orientation[i](1);
158 3 : value.normal[i] = value.normal[i] / value.normal[i].lpNorm<2>();
159 : }
160 : }
161 :
162 6 : amrex::Vector<Set::Scalar> a_eps;
163 : // Diffuse thickness
164 3 : pp_queryarr("eps", a_eps);
165 :
166 3 : if (a_eps.size() < 1)
167 0 : for (int i = 0; i < value.number_of_inclusions; i++)
168 0 : value.eps.push_back(1.e-5);
169 3 : if (a_eps.size() == 1)
170 : {
171 3 : if (a_eps[0] < 0.0)
172 : {
173 0 : Util::Warning(INFO, "Invalid value of laminate.eps. Resetting to 1e-5");
174 0 : a_eps[0] = 1.e-5;
175 : }
176 6 : for (int i = 0; i < value.number_of_inclusions; i++)
177 3 : value.eps.push_back(a_eps[0]);
178 : }
179 : else
180 : {
181 0 : for (int i = 0; i < value.number_of_inclusions; i++)
182 0 : value.eps.push_back(a_eps[i]);
183 : }
184 :
185 6 : std::string mollifier;
186 3 : pp_query("mollifier", mollifier); // Type of mollifer to use (options: dirac, [gaussian])
187 3 : if (mollifier == "Dirac" || mollifier == "dirac")
188 0 : value.moll = Mollifier::Dirac;
189 : else
190 3 : value.moll = Mollifier::Gaussian;
191 :
192 : // Switch to mode where only one component is used.
193 3 : pp_query("singlefab", value.singlefab);
194 :
195 : // Take the complement of the laminate
196 3 : pp_query("invert", value.invert);
197 3 : }
198 : };
199 : }
200 : #endif
|