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