Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/node.h"
22 : #include "libmesh/elem.h"
23 : #include "libmesh/reference_elem.h"
24 : #include "libmesh/libmesh_singleton.h"
25 : #include "libmesh/threads.h"
26 : #include "libmesh/enum_to_string.h"
27 : #include "libmesh/enum_elem_type.h"
28 :
29 : // C++ includes
30 : #include <map>
31 : #include <sstream>
32 : #include <memory> // std::unique_ptr
33 :
34 :
35 : //-----------------------------------------------
36 : // anonymous namespace for implementation details
37 : namespace
38 : {
39 : using namespace libMesh;
40 :
41 : namespace ElemDataStrings
42 : {
43 : // GCC 5.2.0 warns about overlength strings in the auto-generated
44 : // reference_elem.data file.
45 : #pragma GCC diagnostic ignored "-Woverlength-strings"
46 : #include "reference_elem.data"
47 : #pragma GCC diagnostic warning "-Woverlength-strings"
48 : }
49 :
50 : typedef Threads::spin_mutex InitMutex;
51 :
52 : // Mutex for thread safety.
53 : InitMutex init_mtx;
54 :
55 : // map from ElemType to reference element file system object name
56 : typedef std::map<ElemType, const char *> FileMapType;
57 : FileMapType ref_elem_file;
58 : Elem * ref_elem_map[INVALID_ELEM];
59 :
60 :
61 :
62 14 : class SingletonCache : public libMesh::Singleton
63 : {
64 : public:
65 712 : virtual ~SingletonCache() = default;
66 :
67 : std::vector<std::unique_ptr<Node>> node_list;
68 : std::vector<std::unique_ptr<Elem>> elem_list;
69 : };
70 :
71 : // From [0], regarding the lifetime of the singleton_cache variable:
72 : //
73 : // "All variables at namespace level, including the anonymous
74 : // namespace and function local static variable have static
75 : // storage duration unless they are declared thread_local."
76 : //
77 : // Variables with static storage duration are destroyed at the end of
78 : // the program execution. From [1],
79 : //
80 : // "If it is a pointer to the data which is static ... then like all
81 : // other dynamically allocated data, it will only be destructed when
82 : // you delete it. There are two frequent solutions:
83 : // * use a smart pointer, which has a destructor which deletes it, or
84 : // * don't delete it; in most cases, there's really no reason to call the
85 : // destructor, and if you happen to use the instance in the destructors
86 : // of other static objects, you'll run into an order of destruction
87 : // problem."
88 : //
89 : // I tried making the singleton_cache a std::unique_ptr, but this resulted
90 : // in a segfault during program shutdown which appeared to come from the
91 : // single_cache unique_ptr's destructor. I didn't investigate whether the
92 : // issue was caused by a double deletion or what, but it appears that the
93 : // first suggestion above may not be valid in general.
94 : //
95 : // [0]: https://stackoverflow.com/questions/24342393/how-anonymous-namespaces-avoids-making-global-static-variable
96 : // [1]: https://stackoverflow.com/questions/6850009/c-deleting-static-data
97 : SingletonCache * singleton_cache = nullptr;
98 :
99 :
100 :
101 8712 : void read_ref_elem (const ElemType type_in,
102 : std::istream & in)
103 : {
104 336 : libmesh_assert (singleton_cache != nullptr);
105 :
106 336 : std::string dummy;
107 : unsigned int n_elem, n_nodes, elem_type_read, nn;
108 : double x, y, z;
109 :
110 8712 : in >> dummy;
111 8712 : in >> n_elem; /**/ std::getline (in, dummy); libmesh_assert_equal_to (n_elem, 1);
112 8712 : in >> n_nodes; /**/ std::getline (in, dummy);
113 8712 : in >> dummy; /**/ std::getline (in, dummy);
114 8712 : in >> dummy; /**/ std::getline (in, dummy);
115 8712 : in >> dummy; /**/ std::getline (in, dummy);
116 8712 : in >> dummy; /**/ std::getline (in, dummy);
117 8712 : in >> n_elem; /**/ std::getline (in, dummy); libmesh_assert_equal_to (n_elem, 1);
118 :
119 336 : in >> elem_type_read;
120 :
121 336 : libmesh_assert_less (elem_type_read, INVALID_ELEM);
122 336 : libmesh_assert_equal_to (elem_type_read, static_cast<unsigned int>(type_in));
123 336 : libmesh_assert_equal_to (n_nodes, Elem::type_to_n_nodes_map[elem_type_read]);
124 :
125 : // Construct elem of appropriate type, store in the elem_list
126 9048 : auto & uelem = singleton_cache->elem_list.emplace_back(Elem::build(type_in));
127 :
128 : // We are expecting an identity map, so assert it!
129 102729 : for (unsigned int n=0; n<n_nodes; n++)
130 : {
131 3626 : in >> nn;
132 3626 : libmesh_assert_equal_to (n,nn);
133 : }
134 :
135 102729 : for (unsigned int n=0; n<n_nodes; n++)
136 : {
137 3626 : in >> x >> y >> z;
138 :
139 : auto & new_node =
140 97643 : singleton_cache->node_list.emplace_back(Node::build(x,y,z,n));
141 :
142 97643 : uelem->set_node(n, new_node.get());
143 : }
144 :
145 : // it is entirely possible we ran out of file or encountered
146 : // another error. If so, throw an error.
147 9048 : libmesh_error_msg_if(!in, "ERROR while creating element singleton!");
148 :
149 : // Also store a pointer to the newly created Elem in the ref_elem_map array.
150 8376 : ref_elem_map[type_in] = uelem.get();
151 8712 : }
152 :
153 :
154 :
155 134059 : void init_ref_elem_table()
156 : {
157 : // outside mutex - if this pointer is set, we can trust it.
158 134059 : if (singleton_cache != nullptr)
159 11216 : return;
160 :
161 : // playing with fire here - lock before touching shared
162 : // data structures
163 14 : InitMutex::scoped_lock lock(init_mtx);
164 :
165 : // inside mutex - pointer may have changed while waiting
166 : // for the lock to acquire, check it again.
167 363 : if (singleton_cache != nullptr)
168 0 : return;
169 :
170 : // OK, if we get here we have the lock and we are not
171 : // initialized. populate singleton. Note that we do not
172 : // use a smart pointer to manage the singleton_cache variable
173 : // since variables with static storage duration are destroyed
174 : // automatically at the end of program execution.
175 377 : singleton_cache = new SingletonCache;
176 :
177 : // initialize the reference file table
178 : {
179 14 : ref_elem_file.clear();
180 :
181 : // 1D elements
182 363 : ref_elem_file[EDGE2] = ElemDataStrings::one_edge;
183 363 : ref_elem_file[EDGE3] = ElemDataStrings::one_edge3;
184 363 : ref_elem_file[EDGE4] = ElemDataStrings::one_edge4;
185 :
186 : // 2D elements
187 363 : ref_elem_file[TRI3] = ElemDataStrings::one_tri;
188 363 : ref_elem_file[TRI6] = ElemDataStrings::one_tri6;
189 363 : ref_elem_file[TRI7] = ElemDataStrings::one_tri7;
190 :
191 363 : ref_elem_file[QUAD4] = ElemDataStrings::one_quad;
192 363 : ref_elem_file[QUAD8] = ElemDataStrings::one_quad8;
193 363 : ref_elem_file[QUAD9] = ElemDataStrings::one_quad9;
194 :
195 : // 3D elements
196 363 : ref_elem_file[HEX8] = ElemDataStrings::one_hex;
197 363 : ref_elem_file[HEX20] = ElemDataStrings::one_hex20;
198 363 : ref_elem_file[HEX27] = ElemDataStrings::one_hex27;
199 :
200 363 : ref_elem_file[TET4] = ElemDataStrings::one_tet;
201 363 : ref_elem_file[TET10] = ElemDataStrings::one_tet10;
202 363 : ref_elem_file[TET14] = ElemDataStrings::one_tet14;
203 :
204 363 : ref_elem_file[PRISM6] = ElemDataStrings::one_prism;
205 363 : ref_elem_file[PRISM15] = ElemDataStrings::one_prism15;
206 363 : ref_elem_file[PRISM18] = ElemDataStrings::one_prism18;
207 363 : ref_elem_file[PRISM20] = ElemDataStrings::one_prism20;
208 363 : ref_elem_file[PRISM21] = ElemDataStrings::one_prism21;
209 :
210 363 : ref_elem_file[PYRAMID5] = ElemDataStrings::one_pyramid;
211 363 : ref_elem_file[PYRAMID13] = ElemDataStrings::one_pyramid13;
212 363 : ref_elem_file[PYRAMID14] = ElemDataStrings::one_pyramid14;
213 363 : ref_elem_file[PYRAMID18] = ElemDataStrings::one_pyramid18;
214 :
215 : // No entry here for C0POLYGON or C0POLYHEDRON - the reference
216 : // element there depends on the number of nodes and can't be
217 : // precomputed.
218 : }
219 :
220 : // Read'em
221 9075 : for (const auto & [elem_type, filename] : ref_elem_file)
222 : {
223 9720 : std::istringstream stream(filename);
224 8712 : read_ref_elem(elem_type, stream);
225 8040 : }
226 : }
227 :
228 :
229 : // no reason to do this at startup -
230 : // data structures will get initialized *if*
231 : // ReferenceElem::get() is ever called.
232 : // // Class to setup singleton data
233 : // class ReferenceElemSetup : public Singleton::Setup
234 : // {
235 : // void setup ()
236 : // {
237 : // init_ref_elem_table();
238 : // }
239 : // } reference_elem_setup;
240 :
241 : } // anonymous namespace
242 :
243 :
244 :
245 : //----------------------------------------------------------------------------
246 : // external API Implementation
247 : namespace libMesh
248 : {
249 : namespace ReferenceElem
250 : {
251 134059 : const Elem & get (const ElemType type_in)
252 : {
253 11230 : ElemType base_type = type_in;
254 :
255 : // For shell elements, use non shell type as the base type
256 134059 : if (type_in == TRISHELL3)
257 0 : base_type = TRI3;
258 :
259 134059 : if (type_in == QUADSHELL4)
260 3328 : base_type = QUAD4;
261 :
262 100779 : if (type_in == QUADSHELL8)
263 6656 : base_type = QUAD8;
264 :
265 67499 : if (type_in == QUADSHELL9)
266 0 : base_type = QUAD9;
267 :
268 134059 : init_ref_elem_table();
269 :
270 : // Throw an error if the user asked for an ElemType that we don't
271 : // have a reference element for.
272 134059 : libmesh_error_msg_if(ref_elem_map[base_type] == nullptr || type_in == INVALID_ELEM,
273 : "No reference elem data available for ElemType " << type_in
274 : << " = " << Utility::enum_to_string(type_in) << ".");
275 :
276 134059 : return *ref_elem_map[base_type];
277 : }
278 : } // namespace ReferenceElem
279 : } // namespace libMesh
|
