Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #ifndef PALACE_MODELS_MATERIAL_OPERATOR_HPP
5 : #define PALACE_MODELS_MATERIAL_OPERATOR_HPP
6 :
7 : #include <mfem.hpp>
8 : #include "fem/mesh.hpp"
9 : #include "utils/configfile.hpp"
10 :
11 : namespace palace
12 : {
13 :
14 : class IoData;
15 :
16 : //
17 : // A class handling material attributes.
18 : //
19 : class MaterialOperator
20 : {
21 : private:
22 : // Reference to underlying mesh object (not owned).
23 : const Mesh &mesh;
24 :
25 : // Mapping from the local libCEED attribute to material index.
26 : mfem::Array<int> attr_mat;
27 :
28 : // Material properties: relative permeability, relative permittivity, and others (like
29 : // electrical conductivity, London penetration depth for superconductors and Floquet wave
30 : // vector).
31 : mfem::DenseTensor mat_muinv, mat_epsilon, mat_epsilon_imag, mat_epsilon_abs, mat_invz0,
32 : mat_c0, mat_sigma, mat_invLondon, mat_kxTmuinv, mat_muinvkx, mat_kxTmuinvkx, mat_kx;
33 : mfem::DenseMatrix wave_vector_cross;
34 : mfem::Array<double> mat_c0_min, mat_c0_max;
35 :
36 : // Are materials isotropic? True when all the material properties are effectively
37 : // scalar-valued (ie, true scalars or vectors with identical entries). Also true when a
38 : // material is isotropic, the intersection is true when all are isotropic.
39 : mfem::Array<bool> attr_is_isotropic;
40 :
41 : // Flag for global domain attributes with nonzero loss tangent, electrical conductivity,
42 : // London penetration depth, or Floquet wave vector.
43 : bool has_losstan_attr, has_conductivity_attr, has_london_attr, has_wave_attr;
44 :
45 : void SetUpMaterialProperties(const IoData &iodata, const mfem::ParMesh &mesh);
46 : void SetUpFloquetWaveVector(const IoData &iodata, const mfem::ParMesh &mesh);
47 :
48 : // Map from an attribute (specified on a mesh) to a material index (location in the
49 : // property vector).
50 : auto AttrToMat(int attr) const
51 : {
52 152829 : const auto &loc_attr = mesh.GetCeedAttributes();
53 : MFEM_ASSERT(loc_attr.find(attr) != loc_attr.end(),
54 : "Missing libCEED domain attribute for attribute " << attr << "!");
55 124026 : return attr_mat[loc_attr.at(attr) - 1];
56 : }
57 :
58 124020 : auto Wrap(const mfem::DenseTensor &data, int attr) const
59 : {
60 124020 : const int k = AttrToMat(attr);
61 : return mfem::DenseMatrix(const_cast<double *>(data.GetData(k)), data.SizeI(),
62 124020 : data.SizeJ());
63 : }
64 :
65 : public:
66 : MaterialOperator(const IoData &iodata, const Mesh &mesh);
67 :
68 : int SpaceDimension() const { return mat_muinv.SizeI(); }
69 :
70 43092 : auto GetInvPermeability(int attr) const { return Wrap(mat_muinv, attr); }
71 23328 : auto GetPermittivityReal(int attr) const { return Wrap(mat_epsilon, attr); }
72 : auto GetPermittivityImag(int attr) const { return Wrap(mat_epsilon_imag, attr); }
73 : auto GetPermittivityAbs(int attr) const { return Wrap(mat_epsilon_abs, attr); }
74 : auto GetInvImpedance(int attr) const { return Wrap(mat_invz0, attr); }
75 57600 : auto GetLightSpeed(int attr) const { return Wrap(mat_c0, attr); }
76 : auto GetConductivity(int attr) const { return Wrap(mat_sigma, attr); }
77 : auto GetInvLondonDepth(int attr) const { return Wrap(mat_invLondon, attr); }
78 : auto GetFloquetCurl(int attr) const { return Wrap(mat_muinvkx, attr); }
79 : auto GetFloquetMass(int attr) const { return Wrap(mat_kxTmuinvkx, attr); }
80 : auto GetFloquetCross(int attr) const { return Wrap(mat_kx, attr); }
81 :
82 : auto GetLightSpeedMin(int attr) const { return mat_c0_min[AttrToMat(attr)]; }
83 28800 : auto GetLightSpeedMax(int attr) const { return mat_c0_max[AttrToMat(attr)]; }
84 :
85 9 : bool IsIsotropic(int attr) const { return attr_is_isotropic[AttrToMat(attr)]; }
86 :
87 9 : const auto &GetInvPermeability() const { return mat_muinv; }
88 12 : const auto &GetPermittivityReal() const { return mat_epsilon; }
89 0 : const auto &GetPermittivityImag() const { return mat_epsilon_imag; }
90 0 : const auto &GetPermittivityAbs() const { return mat_epsilon_abs; }
91 0 : const auto &GetInvImpedance() const { return mat_invz0; }
92 0 : const auto &GetLightSpeed() const { return mat_c0; }
93 0 : const auto &GetConductivity() const { return mat_sigma; }
94 0 : const auto &GetInvLondonDepth() const { return mat_invLondon; }
95 0 : const auto &GetFloquetCurl() const { return mat_muinvkx; }
96 0 : const auto &GetFloquetMass() const { return mat_kxTmuinvkx; }
97 0 : const auto &GetFloquetCross() const { return mat_kx; }
98 :
99 : const auto &GetLightSpeedMin() const { return mat_c0_min; }
100 0 : const auto &GetLightSpeedMax() const { return mat_c0_max; }
101 :
102 1 : bool HasLossTangent() const { return has_losstan_attr; }
103 0 : bool HasConductivity() const { return has_conductivity_attr; }
104 0 : bool HasLondonDepth() const { return has_london_attr; }
105 0 : bool HasWaveVector() const { return has_wave_attr; }
106 :
107 15 : const auto &GetAttributeToMaterial() const { return attr_mat; }
108 : mfem::Array<int> GetBdrAttributeToMaterial() const;
109 :
110 : template <typename T>
111 : auto GetCeedAttributes(const T &attr_list) const
112 : {
113 0 : return mesh.GetCeedAttributes(attr_list);
114 : }
115 : template <typename T>
116 : auto GetCeedBdrAttributes(const T &attr_list) const
117 : {
118 0 : return mesh.GetCeedBdrAttributes(attr_list);
119 : }
120 :
121 6 : auto MaxCeedAttribute() const { return mesh.MaxCeedAttribute(); }
122 0 : auto MaxCeedBdrAttribute() const { return mesh.MaxCeedBdrAttribute(); }
123 :
124 : const auto &GetMesh() const { return mesh; }
125 : };
126 :
127 : //
128 : // Material property represented as a piecewise constant coefficient over domain or boundary
129 : // mesh elements. Can be scalar-valued or matrix-valued. This should probably always operate
130 : // at the level of libCEED attribute numbers (contiguous, 1-based) for consistency.
131 : //
132 : class MaterialPropertyCoefficient
133 : {
134 : private:
135 : // Map attribute to material index (coeff = mat_coeff[attr_mat[attr - 1]], for 1-based
136 : // attributes).
137 : mfem::Array<int> attr_mat;
138 :
139 : // Material propetry coefficients, ordered by material index.
140 : mfem::DenseTensor mat_coeff;
141 :
142 : public:
143 : MaterialPropertyCoefficient(int attr_max);
144 : MaterialPropertyCoefficient(const mfem::Array<int> &attr_mat_,
145 : const mfem::DenseTensor &mat_coeff_, double a = 1.0);
146 :
147 0 : bool empty() const { return mat_coeff.TotalSize() == 0; }
148 :
149 9934 : const auto &GetAttributeToMaterial() const { return attr_mat; }
150 : const auto &GetMaterialProperties() const { return mat_coeff; }
151 :
152 : void AddCoefficient(const mfem::Array<int> &attr_mat_,
153 : const mfem::DenseTensor &mat_coeff_, double a = 1.0);
154 :
155 : template <typename T>
156 : void AddMaterialProperty(const mfem::Array<int> &attr_list, const T &coeff,
157 : double a = 1.0);
158 : template <typename T>
159 : void AddMaterialProperty(int attr, const T &coeff, double a = 1.0)
160 : {
161 : mfem::Array<int> attr_list(1);
162 : attr_list[0] = attr;
163 : AddMaterialProperty(attr_list, coeff, a);
164 : }
165 :
166 : MaterialPropertyCoefficient &operator*=(double a);
167 :
168 : void RestrictCoefficient(const mfem::Array<int> &attr_list);
169 :
170 : void NormalProjectedCoefficient(const mfem::Vector &normal);
171 : };
172 :
173 : } // namespace palace
174 :
175 : namespace palace::internal::mat
176 : {
177 :
178 : template <std::size_t N>
179 : bool IsOrthonormal(const config::SymmetricMatrixData<N> &data);
180 :
181 : template <std::size_t N>
182 : bool IsValid(const config::SymmetricMatrixData<N> &data);
183 :
184 : template <std::size_t N>
185 : bool IsIsotropic(const config::SymmetricMatrixData<N> &data);
186 :
187 : template <std::size_t N>
188 : bool IsIdentity(const config::SymmetricMatrixData<N> &data);
189 :
190 : } // namespace palace::internal::mat
191 :
192 : #endif // PALACE_MODELS_MATERIAL_OPERATOR_HPP
|
