Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "materialoperator.hpp"
5 :
6 : #include <cmath>
7 : #include <functional>
8 : #include <limits>
9 : #include "linalg/densematrix.hpp"
10 : #include "utils/communication.hpp"
11 : #include "utils/geodata.hpp"
12 : #include "utils/iodata.hpp"
13 :
14 : namespace palace
15 : {
16 :
17 : namespace internal::mat
18 : {
19 :
20 : template <std::size_t N>
21 272 : bool IsOrthonormal(const config::SymmetricMatrixData<N> &data)
22 : {
23 :
24 : // All the vectors are normalized.
25 : constexpr auto tol = 1.0e-6;
26 : auto UnitNorm = [&](const std::array<double, N> &x)
27 : {
28 : double s = -1.0;
29 3244 : for (const auto &i : x)
30 : {
31 2433 : s += std::pow(i, 2);
32 : }
33 : return std::abs(s) < tol;
34 : };
35 : bool valid = std::all_of(data.v.begin(), data.v.end(), UnitNorm);
36 :
37 : // All the vectors are orthogonal.
38 1088 : for (std::size_t i1 = 0; i1 < N; i1++)
39 : {
40 : const auto &v1 = data.v.at(i1);
41 1632 : for (std::size_t i2 = i1 + 1; i2 < N; i2++)
42 : {
43 : const auto &v2 = data.v.at(i2);
44 : double s = 0.0;
45 3264 : for (std::size_t j = 0; j < N; j++)
46 : {
47 2448 : s += v1[j] * v2[j];
48 : }
49 816 : valid &= std::abs(s) < tol;
50 : }
51 : }
52 272 : return valid;
53 : }
54 :
55 : template <std::size_t N>
56 114 : bool IsValid(const config::SymmetricMatrixData<N> &data)
57 : {
58 114 : return IsOrthonormal(data) && std::all_of(data.s.begin(), data.s.end(),
59 114 : [](auto d) { return std::abs(d) > 0.0; });
60 : }
61 :
62 : template <std::size_t N>
63 146 : bool IsIsotropic(const config::SymmetricMatrixData<N> &data)
64 : {
65 146 : return IsOrthonormal(data) &&
66 430 : std::all_of(data.s.begin(), data.s.end(), [&](auto d) { return d == data.s[0]; });
67 : }
68 :
69 : template <std::size_t N>
70 9 : bool IsIdentity(const config::SymmetricMatrixData<N> &data)
71 : {
72 9 : return IsOrthonormal(data) &&
73 9 : std::all_of(data.s.begin(), data.s.end(), [](auto d) { return d == 1.0; });
74 : }
75 :
76 : template <std::size_t N>
77 190 : mfem::DenseMatrix ToDenseMatrix(const config::SymmetricMatrixData<N> &data)
78 : {
79 190 : mfem::DenseMatrix M(N, N);
80 : mfem::Vector V(N);
81 760 : for (std::size_t i = 0; i < N; i++)
82 : {
83 2280 : for (std::size_t j = 0; j < N; j++)
84 : {
85 1710 : V(j) = data.v[i][j];
86 : }
87 570 : AddMult_a_VVt(data.s[i], V, M);
88 : }
89 190 : return M;
90 0 : }
91 :
92 : } // namespace internal::mat
93 :
94 78 : MaterialOperator::MaterialOperator(const IoData &iodata, const Mesh &mesh) : mesh(mesh)
95 : {
96 39 : SetUpMaterialProperties(iodata, mesh);
97 39 : }
98 :
99 39 : void MaterialOperator::SetUpMaterialProperties(const IoData &iodata,
100 : const mfem::ParMesh &mesh)
101 : {
102 : // Check that material attributes have been specified correctly. The mesh attributes may
103 : // be non-contiguous and when no material attribute is specified the elements are deleted
104 : // from the mesh so as to not cause problems.
105 39 : MFEM_VERIFY(!iodata.domains.materials.empty(), "Materials must be non-empty!");
106 : {
107 39 : int attr_max = mesh.attributes.Size() ? mesh.attributes.Max() : 0;
108 : mfem::Array<int> attr_marker(attr_max);
109 : attr_marker = 0;
110 78 : for (auto attr : mesh.attributes)
111 : {
112 39 : attr_marker[attr - 1] = 1;
113 : }
114 78 : for (const auto &data : iodata.domains.materials)
115 : {
116 77 : for (auto attr : data.attributes)
117 : {
118 38 : MFEM_VERIFY(
119 : attr > 0 && attr <= attr_max,
120 : "Material attribute tags must be non-negative and correspond to attributes "
121 : "in the mesh!");
122 38 : MFEM_VERIFY(attr_marker[attr - 1], "Unknown material attribute " << attr << "!");
123 : }
124 : }
125 : }
126 :
127 : // Set up material properties of the different domain regions, represented with element-
128 : // wise constant matrix-valued coefficients for the relative permeability, permittivity,
129 : // and other material properties.
130 39 : const auto &loc_attr = this->mesh.GetCeedAttributes();
131 39 : mfem::Array<int> mat_marker(iodata.domains.materials.size());
132 : mat_marker = 0;
133 : int nmats = 0;
134 78 : for (std::size_t i = 0; i < iodata.domains.materials.size(); i++)
135 : {
136 39 : const auto &data = iodata.domains.materials[i];
137 39 : for (auto attr : data.attributes)
138 : {
139 38 : if (loc_attr.find(attr) != loc_attr.end())
140 : {
141 38 : mat_marker[i] = 1;
142 38 : nmats++;
143 38 : break;
144 : }
145 : }
146 : }
147 39 : attr_mat.SetSize(loc_attr.size());
148 : attr_mat = -1;
149 :
150 39 : attr_is_isotropic.SetSize(nmats);
151 :
152 : const int sdim = mesh.SpaceDimension();
153 39 : mat_muinv.SetSize(sdim, sdim, nmats);
154 39 : mat_epsilon.SetSize(sdim, sdim, nmats);
155 39 : mat_epsilon_imag.SetSize(sdim, sdim, nmats);
156 39 : mat_epsilon_abs.SetSize(sdim, sdim, nmats);
157 39 : mat_invz0.SetSize(sdim, sdim, nmats);
158 39 : mat_c0.SetSize(sdim, sdim, nmats);
159 39 : mat_sigma.SetSize(sdim, sdim, nmats);
160 39 : mat_invLondon.SetSize(sdim, sdim, nmats);
161 39 : mat_c0_min.SetSize(nmats);
162 39 : mat_c0_max.SetSize(nmats);
163 39 : mat_muinvkx.SetSize(sdim, sdim, nmats);
164 39 : mat_kxTmuinvkx.SetSize(sdim, sdim, nmats);
165 39 : mat_kx.SetSize(sdim, sdim, nmats);
166 39 : has_losstan_attr = has_conductivity_attr = has_london_attr = has_wave_attr = false;
167 :
168 : // Set up Floquet wave vector for periodic meshes with phase-delay constraints.
169 39 : SetUpFloquetWaveVector(iodata, mesh);
170 :
171 : int count = 0;
172 78 : for (std::size_t i = 0; i < iodata.domains.materials.size(); i++)
173 : {
174 39 : if (!mat_marker[i])
175 : {
176 1 : continue;
177 : }
178 : const auto &data = iodata.domains.materials[i];
179 38 : if (iodata.problem.type == ProblemType::ELECTROSTATIC)
180 : {
181 3 : MFEM_VERIFY(internal::mat::IsValid(data.epsilon_r),
182 : "Material has no valid permittivity defined!");
183 3 : if (!internal::mat::IsIdentity(data.mu_r) || internal::mat::IsValid(data.sigma) ||
184 3 : std::abs(data.lambda_L) > 0.0)
185 : {
186 0 : Mpi::Warning(
187 : "Electrostatic problem type does not account for material permeability,\n"
188 : "electrical conductivity, or London depth!\n");
189 : }
190 : }
191 35 : else if (iodata.problem.type == ProblemType::MAGNETOSTATIC)
192 : {
193 3 : MFEM_VERIFY(internal::mat::IsValid(data.mu_r),
194 : "Material has no valid permeability defined!");
195 6 : if (!internal::mat::IsIdentity(data.epsilon_r) ||
196 6 : internal::mat::IsValid(data.tandelta) || internal::mat::IsValid(data.sigma) ||
197 3 : std::abs(data.lambda_L) > 0.0)
198 : {
199 0 : Mpi::Warning(
200 : "Magnetostatic problem type does not account for material permittivity,\n"
201 : "loss tangent, electrical conductivity, or London depth!\n");
202 : }
203 : }
204 : else
205 : {
206 32 : MFEM_VERIFY(internal::mat::IsValid(data.mu_r) &&
207 : internal::mat::IsValid(data.epsilon_r),
208 : "Material has no valid permeability or no valid permittivity defined!");
209 32 : if (iodata.problem.type == ProblemType::TRANSIENT)
210 : {
211 3 : MFEM_VERIFY(!internal::mat::IsValid(data.tandelta),
212 : "Transient problem type does not support material loss tangent, use "
213 : "electrical conductivity instead!");
214 : }
215 : else
216 : {
217 29 : MFEM_VERIFY(
218 : !(internal::mat::IsValid(data.tandelta) && internal::mat::IsValid(data.sigma)),
219 : "Material loss model should probably use only one of loss tangent or "
220 : "electrical conductivity!");
221 : }
222 : }
223 :
224 112 : attr_is_isotropic[i] = internal::mat::IsIsotropic(data.mu_r) &&
225 71 : internal::mat::IsIsotropic(data.epsilon_r) &&
226 107 : internal::mat::IsIsotropic(data.tandelta) &&
227 34 : internal::mat::IsIsotropic(data.sigma);
228 :
229 : // Map all attributes to this material property index.
230 76 : for (auto attr : data.attributes)
231 : {
232 : auto it = loc_attr.find(attr);
233 38 : if (it != loc_attr.end())
234 : {
235 38 : MFEM_VERIFY(
236 : attr_mat[it->second - 1] < 0,
237 : "Detected multiple definitions of material properties for domain attribute "
238 : << attr << "!");
239 38 : attr_mat[it->second - 1] = count;
240 : }
241 : }
242 :
243 : // Compute the inverse of the input permeability matrix.
244 38 : mfem::DenseMatrix mat_mu = internal::mat::ToDenseMatrix(data.mu_r);
245 76 : mfem::DenseMatrixInverse(mat_mu, true).GetInverseMatrix(mat_muinv(count));
246 :
247 : // Material permittivity: Re{ε} = ε, Im{ε} = -ε * tan(δ)
248 38 : mfem::DenseMatrix T(sdim, sdim);
249 76 : mat_epsilon(count) = internal::mat::ToDenseMatrix(data.epsilon_r);
250 76 : Mult(mat_epsilon(count), internal::mat::ToDenseMatrix(data.tandelta), T);
251 38 : T *= -1.0;
252 38 : mat_epsilon_imag(count) = T;
253 38 : if (mat_epsilon_imag(count).MaxMaxNorm() > 0.0)
254 : {
255 1 : has_losstan_attr = true;
256 : }
257 :
258 : // ε * √(I + tan(δ) * tan(δ)ᵀ)
259 38 : MultAAt(internal::mat::ToDenseMatrix(data.tandelta), T);
260 152 : for (int d = 0; d < T.Height(); d++)
261 : {
262 114 : T(d, d) += 1.0;
263 : }
264 76 : Mult(mat_epsilon(count), linalg::MatrixSqrt(T), mat_epsilon_abs(count));
265 :
266 : // √(μ⁻¹ ε)
267 38 : Mult(mat_muinv(count), mat_epsilon(count), mat_invz0(count));
268 76 : mat_invz0(count) = linalg::MatrixSqrt(mat_invz0(count));
269 :
270 : // √((μ ε)⁻¹)
271 38 : Mult(mat_mu, mat_epsilon(count), T);
272 76 : mat_c0(count) = linalg::MatrixPow(T, -0.5);
273 38 : mat_c0_min[count] = linalg::SingularValueMin(mat_c0(count));
274 38 : mat_c0_max[count] = linalg::SingularValueMax(mat_c0(count));
275 :
276 : // Electrical conductivity, σ
277 76 : mat_sigma(count) = internal::mat::ToDenseMatrix(data.sigma);
278 38 : if (mat_sigma(count).MaxMaxNorm() > 0.0)
279 : {
280 1 : has_conductivity_attr = true;
281 : }
282 :
283 : // λ⁻² * μ⁻¹
284 38 : mat_invLondon(count) = mat_muinv(count);
285 : mat_invLondon(count) *=
286 76 : std::abs(data.lambda_L) > 0.0 ? std::pow(data.lambda_L, -2.0) : 0.0;
287 38 : if (mat_invLondon(count).MaxMaxNorm() > 0.0)
288 : {
289 0 : has_london_attr = true;
290 : }
291 :
292 : // μ⁻¹ [k x]
293 38 : Mult(mat_muinv(count), wave_vector_cross, mat_muinvkx(count));
294 :
295 : // [k x]^T μ⁻¹ [k x]
296 38 : T.Transpose(wave_vector_cross);
297 38 : Mult(T, mat_muinvkx(count), mat_kxTmuinvkx(count));
298 :
299 : // [k x]
300 38 : mat_kx(count) = wave_vector_cross;
301 :
302 38 : count++;
303 38 : }
304 39 : bool has_attr[4] = {has_losstan_attr, has_conductivity_attr, has_london_attr,
305 39 : has_wave_attr};
306 : Mpi::GlobalOr(4, has_attr, mesh.GetComm());
307 39 : has_losstan_attr = has_attr[0];
308 39 : has_conductivity_attr = has_attr[1];
309 39 : has_london_attr = has_attr[2];
310 39 : has_wave_attr = has_attr[3];
311 39 : }
312 :
313 39 : void MaterialOperator::SetUpFloquetWaveVector(const IoData &iodata,
314 : const mfem::ParMesh &mesh)
315 : {
316 : const int sdim = mesh.SpaceDimension();
317 : const double tol = std::numeric_limits<double>::epsilon();
318 :
319 : // Get Floquet wave vector.
320 : mfem::Vector wave_vector(sdim);
321 39 : wave_vector = 0.0;
322 : const auto &data = iodata.boundaries.periodic;
323 39 : MFEM_VERIFY(static_cast<int>(data.wave_vector.size()) == sdim,
324 : "Floquet wave vector size must equal the spatial dimension.");
325 : std::copy(data.wave_vector.begin(), data.wave_vector.end(), wave_vector.GetData());
326 39 : has_wave_attr = (wave_vector.Norml2() > tol);
327 :
328 39 : MFEM_VERIFY(!has_wave_attr || iodata.problem.type == ProblemType::DRIVEN ||
329 : iodata.problem.type == ProblemType::EIGENMODE,
330 : "Quasi-periodic Floquet boundary conditions are only available for "
331 : " frequency domain driven or eigenmode simulations!");
332 : MFEM_VERIFY(!has_wave_attr || sdim == 3,
333 : "Quasi-periodic Floquet periodic boundary conditions are only available "
334 : " in 3D!");
335 :
336 : // Get mesh dimensions in x/y/z coordinates.
337 : mfem::Vector bbmin, bbmax;
338 39 : mesh::GetAxisAlignedBoundingBox(mesh, bbmin, bbmax);
339 39 : bbmax -= bbmin;
340 :
341 : // Ensure Floquet wave vector components are in range [-π/L, π/L].
342 156 : for (int i = 0; i < sdim; i++)
343 : {
344 117 : if (wave_vector[i] > M_PI / bbmax[i])
345 : {
346 0 : wave_vector[i] =
347 0 : -M_PI / bbmax[i] + fmod(wave_vector[i] + M_PI / bbmax[i], 2 * M_PI / bbmax[i]);
348 : }
349 117 : else if (wave_vector[i] < M_PI / bbmax[i])
350 : {
351 117 : wave_vector[i] =
352 117 : M_PI / bbmax[i] + fmod(wave_vector[i] - M_PI / bbmax[i], 2 * M_PI / bbmax[i]);
353 : }
354 : }
355 :
356 : // Matrix representation of cross product with wave vector
357 : // [k x] = | 0 -k3 k2|
358 : // | k3 0 -k1|
359 : // |-k2 k1 0 |
360 39 : wave_vector_cross.SetSize(3);
361 39 : wave_vector_cross = 0.0;
362 39 : wave_vector_cross(0, 1) = -wave_vector[2];
363 39 : wave_vector_cross(0, 2) = wave_vector[1];
364 39 : wave_vector_cross(1, 0) = wave_vector[2];
365 39 : wave_vector_cross(1, 2) = -wave_vector[0];
366 39 : wave_vector_cross(2, 0) = -wave_vector[1];
367 39 : wave_vector_cross(2, 1) = wave_vector[0];
368 39 : }
369 :
370 0 : mfem::Array<int> MaterialOperator::GetBdrAttributeToMaterial() const
371 : {
372 : // Construct map from all (contiguous) local libCEED boundary attributes to the material
373 : // index in the neighboring element.
374 0 : mfem::Array<int> bdr_attr_mat(mesh.MaxCeedBdrAttribute());
375 : bdr_attr_mat = -1;
376 0 : for (const auto &[attr, bdr_attr_map] : mesh.GetCeedBdrAttributes())
377 : {
378 0 : for (auto it = bdr_attr_map.begin(); it != bdr_attr_map.end(); ++it)
379 : {
380 : MFEM_ASSERT(it->second > 0 && it->second <= bdr_attr_mat.Size(),
381 : "Invalid libCEED boundary attribute " << it->second << "!");
382 0 : bdr_attr_mat[it->second - 1] = AttrToMat(it->first);
383 : }
384 : }
385 0 : return bdr_attr_mat;
386 : }
387 :
388 3684 : MaterialPropertyCoefficient::MaterialPropertyCoefficient(int attr_max)
389 : {
390 3684 : attr_mat.SetSize(attr_max);
391 : attr_mat = -1;
392 3684 : }
393 :
394 7432 : MaterialPropertyCoefficient::MaterialPropertyCoefficient(
395 7432 : const mfem::Array<int> &attr_mat_, const mfem::DenseTensor &mat_coeff_, double a)
396 7432 : : attr_mat(attr_mat_), mat_coeff(mat_coeff_)
397 : {
398 7432 : *this *= a;
399 7432 : }
400 :
401 : namespace
402 : {
403 :
404 0 : void UpdateProperty(mfem::DenseTensor &mat_coeff, int k, double coeff, double a)
405 : {
406 : // Constant diagonal coefficient.
407 0 : if (mat_coeff.SizeI() == 0 && mat_coeff.SizeJ() == 0)
408 : {
409 : // Initialize the coefficient material properties.
410 0 : MFEM_VERIFY(k == 0 && mat_coeff.SizeK() == 1,
411 : "Unexpected initial size for MaterialPropertyCoefficient!");
412 0 : mat_coeff.SetSize(1, 1, mat_coeff.SizeK());
413 0 : mat_coeff(0, 0, k) = a * coeff;
414 : }
415 : else
416 : {
417 0 : MFEM_VERIFY(mat_coeff.SizeI() == mat_coeff.SizeJ(),
418 : "Invalid dimensions for MaterialPropertyCoefficient update!");
419 0 : for (int i = 0; i < mat_coeff.SizeI(); i++)
420 : {
421 0 : mat_coeff(i, i, k) += a * coeff;
422 : }
423 : }
424 0 : }
425 :
426 0 : void UpdateProperty(mfem::DenseTensor &mat_coeff, int k, const mfem::DenseMatrix &coeff,
427 : double a)
428 : {
429 0 : if (mat_coeff.SizeI() == 0 && mat_coeff.SizeJ() == 0)
430 : {
431 : // Initialize the coefficient material properties.
432 0 : MFEM_VERIFY(k == 0 && mat_coeff.SizeK() == 1,
433 : "Unexpected initial size for MaterialPropertyCoefficient!");
434 0 : mat_coeff.SetSize(coeff.Height(), coeff.Width(), mat_coeff.SizeK());
435 0 : mat_coeff(k).Set(a, coeff);
436 : }
437 0 : else if (coeff.Height() == mat_coeff.SizeI() && coeff.Width() == mat_coeff.SizeJ())
438 : {
439 : // Add as full matrix.
440 0 : mat_coeff(k).Add(a, coeff);
441 : }
442 0 : else if (coeff.Height() == 1 && coeff.Width() == 1)
443 : {
444 : // Add as diagonal.
445 0 : UpdateProperty(mat_coeff, k, coeff(0, 0), a);
446 : }
447 0 : else if (mat_coeff.SizeI() == 1 && mat_coeff.SizeJ() == 1)
448 : {
449 : // Convert to matrix coefficient and previous data add as diagonal.
450 0 : mfem::DenseTensor mat_coeff_scalar(mat_coeff);
451 0 : mat_coeff.SetSize(coeff.Height(), coeff.Width(), mat_coeff_scalar.SizeK());
452 0 : mat_coeff = 0.0;
453 0 : for (int l = 0; l < mat_coeff.SizeK(); l++)
454 : {
455 0 : UpdateProperty(mat_coeff, l, mat_coeff_scalar(0, 0, l), 1.0);
456 : }
457 0 : mat_coeff(k).Add(a, coeff);
458 : }
459 : else
460 : {
461 0 : MFEM_ABORT("Invalid dimensions when updating material property at index " << k << "!");
462 : }
463 0 : }
464 :
465 0 : bool Equals(const mfem::DenseMatrix &mat_coeff, double coeff, double a)
466 : {
467 0 : MFEM_VERIFY(mat_coeff.Height() == mat_coeff.Width(),
468 : "Invalid dimensions for MaterialPropertyCoefficient update!");
469 : constexpr double tol = 1.0e-9;
470 0 : for (int i = 0; i < mat_coeff.Height(); i++)
471 : {
472 0 : if (std::abs(mat_coeff(i, i) - a * coeff) >= tol * std::abs(mat_coeff(i, i)))
473 : {
474 : return false;
475 : }
476 0 : for (int j = 0; j < mat_coeff.Width(); j++)
477 : {
478 0 : if (j != i && std::abs(mat_coeff(i, j)) > 0.0)
479 : {
480 : return false;
481 : }
482 : }
483 : }
484 : return true;
485 : }
486 :
487 0 : bool Equals(const mfem::DenseMatrix &mat_coeff, const mfem::DenseMatrix &coeff, double a)
488 : {
489 0 : if (coeff.Height() == 1 && coeff.Width() == 1)
490 : {
491 0 : return Equals(mat_coeff, coeff(0, 0), a);
492 : }
493 : else
494 : {
495 : constexpr double tol = 1.0e-9;
496 0 : mfem::DenseMatrix T(mat_coeff);
497 0 : T.Add(-a, coeff);
498 0 : return (T.MaxMaxNorm() < tol * mat_coeff.MaxMaxNorm());
499 0 : }
500 : }
501 :
502 : } // namespace
503 :
504 12 : void MaterialPropertyCoefficient::AddCoefficient(const mfem::Array<int> &attr_mat_,
505 : const mfem::DenseTensor &mat_coeff_,
506 : double a)
507 : {
508 12 : if (empty())
509 : {
510 12 : MFEM_VERIFY(attr_mat_.Size() == attr_mat.Size(),
511 : "Invalid resize of attribute to material property map in "
512 : "MaterialPropertyCoefficient::AddCoefficient!");
513 12 : attr_mat = attr_mat_;
514 12 : mat_coeff = mat_coeff_;
515 12 : *this *= a;
516 : }
517 0 : else if (attr_mat_ == attr_mat)
518 : {
519 0 : MFEM_VERIFY(mat_coeff_.SizeK() == mat_coeff.SizeK(),
520 : "Invalid dimensions for MaterialPropertyCoefficient::AddCoefficient!");
521 0 : for (int k = 0; k < mat_coeff.SizeK(); k++)
522 : {
523 0 : UpdateProperty(mat_coeff, k, mat_coeff_(k), a);
524 : }
525 : }
526 : else
527 : {
528 0 : for (int k = 0; k < mat_coeff_.SizeK(); k++)
529 : {
530 : // Get list of all attributes which use this material property.
531 : mfem::Array<int> attr_list;
532 : attr_list.Reserve(attr_mat_.Size());
533 0 : for (int i = 0; i < attr_mat_.Size(); i++)
534 : {
535 0 : if (attr_mat_[i] == k)
536 : {
537 0 : attr_list.Append(i + 1);
538 : }
539 : }
540 :
541 : // Add or update the material property.
542 0 : AddMaterialProperty(attr_list, mat_coeff_(k), a);
543 : }
544 : }
545 12 : }
546 :
547 : template <typename T>
548 0 : void MaterialPropertyCoefficient::AddMaterialProperty(const mfem::Array<int> &attr_list,
549 : const T &coeff, double a)
550 : {
551 : // Preprocess the attribute list. If any of the given attributes already have material
552 : // properties assigned, then they all need to point to the same material and it is
553 : // updated in place. Otherwise a new material is added for these attributes.
554 0 : if (attr_list.Size() == 0)
555 : {
556 : // No attributes, nothing to add.
557 : return;
558 : }
559 :
560 : int mat_idx = -1;
561 0 : for (auto attr : attr_list)
562 : {
563 0 : MFEM_VERIFY(attr <= attr_mat.Size(),
564 : "Out of bounds access for attribute "
565 : << attr << " in MaterialPropertyCoefficient::AddMaterialProperty!");
566 0 : if (mat_idx < 0)
567 : {
568 0 : mat_idx = attr_mat[attr - 1];
569 : }
570 : else
571 : {
572 0 : MFEM_VERIFY(mat_idx == attr_mat[attr - 1],
573 : "All attributes for MaterialPropertyCoefficient::AddMaterialProperty "
574 : "must correspond to the same "
575 : "existing material if it exists!");
576 : }
577 : }
578 :
579 0 : if (mat_idx < 0)
580 : {
581 : // Check if we can reuse an existing material.
582 0 : for (int k = 0; k < mat_coeff.SizeK(); k++)
583 : {
584 0 : if (Equals(mat_coeff(k), coeff, a))
585 : {
586 : mat_idx = k;
587 : break;
588 : }
589 : }
590 0 : if (mat_idx < 0)
591 : {
592 : // Append a new material and assign the attributes to it.
593 0 : const mfem::DenseTensor mat_coeff_backup(mat_coeff);
594 0 : mat_coeff.SetSize(mat_coeff_backup.SizeI(), mat_coeff_backup.SizeJ(),
595 : mat_coeff_backup.SizeK() + 1);
596 0 : for (int k = 0; k < mat_coeff_backup.SizeK(); k++)
597 : {
598 0 : mat_coeff(k) = mat_coeff_backup(k);
599 : }
600 0 : mat_idx = mat_coeff.SizeK() - 1;
601 : }
602 0 : mat_coeff(mat_idx) = 0.0; // Zero out so we can add
603 :
604 : // Assign all attributes to this new material.
605 0 : for (auto attr : attr_list)
606 : {
607 0 : attr_mat[attr - 1] = mat_idx;
608 : }
609 : }
610 0 : UpdateProperty(mat_coeff, mat_idx, coeff, a);
611 : }
612 :
613 7444 : MaterialPropertyCoefficient &MaterialPropertyCoefficient::operator*=(double a)
614 : {
615 37112 : for (int k = 0; k < mat_coeff.SizeK(); k++)
616 : {
617 29668 : mat_coeff(k) *= a;
618 : }
619 7444 : return *this;
620 : }
621 :
622 0 : void MaterialPropertyCoefficient::RestrictCoefficient(const mfem::Array<int> &attr_list)
623 : {
624 : // Create a new material property coefficient with materials corresponding to only the
625 : // unique ones in the given attribute list.
626 0 : const mfem::Array<int> attr_mat_orig(attr_mat);
627 0 : const mfem::DenseTensor mat_coeff_orig(mat_coeff);
628 : attr_mat = -1;
629 0 : mat_coeff.SetSize(mat_coeff_orig.SizeI(), mat_coeff_orig.SizeJ(), 0);
630 0 : for (auto attr : attr_list)
631 : {
632 0 : if (attr_mat[attr - 1] >= 0)
633 : {
634 : // Attribute has already been processed.
635 0 : continue;
636 : }
637 :
638 : // Find all attributes in restricted list of attributes which map to this material index
639 : // and process them together.
640 0 : const int orig_mat_idx = attr_mat_orig[attr - 1];
641 : const int new_mat_idx = mat_coeff.SizeK();
642 0 : for (auto attr2 : attr_list)
643 : {
644 0 : if (attr_mat_orig[attr2 - 1] == orig_mat_idx)
645 : {
646 0 : attr_mat[attr2 - 1] = new_mat_idx;
647 : }
648 : }
649 :
650 : // Append the new material property.
651 0 : const mfem::DenseTensor mat_coeff_backup(mat_coeff);
652 0 : mat_coeff.SetSize(mat_coeff_backup.SizeI(), mat_coeff_backup.SizeJ(),
653 : mat_coeff_backup.SizeK() + 1);
654 0 : for (int k = 0; k < mat_coeff_backup.SizeK(); k++)
655 : {
656 0 : mat_coeff(k) = mat_coeff_backup(k);
657 : }
658 0 : mat_coeff(new_mat_idx) = mat_coeff_orig(orig_mat_idx);
659 : }
660 0 : }
661 :
662 0 : void MaterialPropertyCoefficient::NormalProjectedCoefficient(const mfem::Vector &normal)
663 : {
664 0 : mfem::DenseTensor mat_coeff_backup(mat_coeff);
665 0 : mat_coeff.SetSize(1, 1, mat_coeff_backup.SizeK());
666 0 : for (int k = 0; k < mat_coeff.SizeK(); k++)
667 : {
668 0 : mat_coeff(k) = mat_coeff_backup(k).InnerProduct(normal, normal);
669 : }
670 0 : }
671 :
672 : template void MaterialPropertyCoefficient::AddMaterialProperty(const mfem::Array<int> &,
673 : const mfem::DenseMatrix &,
674 : double);
675 : template void MaterialPropertyCoefficient::AddMaterialProperty(const mfem::Array<int> &,
676 : const double &, double);
677 :
678 : // Explicit template instantiations for internal::mat functions.
679 : template bool internal::mat::IsOrthonormal(const config::SymmetricMatrixData<3> &);
680 : template bool internal::mat::IsValid(const config::SymmetricMatrixData<3> &);
681 : template bool internal::mat::IsIsotropic(const config::SymmetricMatrixData<3> &);
682 : template bool internal::mat::IsIdentity(const config::SymmetricMatrixData<3> &);
683 :
684 : } // namespace palace
|
