1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "mfem.hpp"
13 #include "unit_tests.hpp"
14 #include <fstream>
15 #include <sstream>
16
17 using namespace mfem;
18
19 /**
20 * \file test_linear_fes
21 *
22 * Some simple tests to confirm that linear finite element collections
23 * and H1_FECollections of order 1 have the correct number of dofs.
24 * Also optionally outputs the mesh and grid functions in mfem and vtk formats
25 */
26
27 namespace
28 {
29 static const bool dumpMeshes = true; // will dump mesh files when true
30 }
31
32 /// Generate a name for the file
33 /// using the mesh's element type, the fec collection and the file extension
34 template<typename FEColType>
generateMeshFilename(FEColType & fec,Mesh & mesh,const std::string & prefix="",const std::string & fileExt="")35 std::string generateMeshFilename(FEColType& fec, Mesh& mesh,
36 const std::string& prefix = "",
37 const std::string& fileExt = "")
38 {
39 std::stringstream sstr;
40
41 REQUIRE ( mesh.GetNE() == 1 );
42 int geom = mesh.GetElementBaseGeometry(0);
43
44 sstr << prefix
45 << Geometry::Name[geom]<< "_" << fec.Name()
46 << fileExt;
47
48 return sstr.str();
49 }
50
51
52
53 /// Tests a scalar and vector grid function on the fes defined by fec and mesh
54 template<typename FEColType>
testGridFunctions(FEColType & fec,Mesh & mesh,int expScalarDofs)55 void testGridFunctions(FEColType& fec, Mesh& mesh, int expScalarDofs)
56 {
57 // Note: we assume the mesh has a single element
58 REQUIRE( mesh.GetNE() == 1);
59 const int eltId = 0;
60 const int geom = mesh.GetElementBaseGeometry(eltId);
61
62 // Create a scalar grid function on this fec and test setting/getting values
63 FiniteElementSpace scalarFes(&mesh, &fec);
64 GridFunction scalarGF(&scalarFes);
65 {
66 GridFunction& gf = scalarGF;
67 REQUIRE( gf.FESpace()->GetVDim() == 1 );
68 REQUIRE( gf.Size() == expScalarDofs );
69 gf = 0.;
70
71 // Set the values of the dofs to non-zero values
72 mfem::Array<int> dofs;
73
74 gf.FESpace()->GetElementDofs(eltId,dofs);
75 REQUIRE( dofs.Size() > 0 );
76
77 for (int i = 0; i < dofs.Size(); ++i)
78 {
79 gf[ dofs[i] ] = (1. + i)/(1. + dofs.Size());
80 }
81
82 // Test access to the dof values by getting the value at the elt midpoint
83 double centerValue = gf.GetValue( eltId, Geometries.GetCenter(geom) );
84 REQUIRE( centerValue > 0. );
85 REQUIRE( centerValue < 1. );
86 }
87
88 // Create a vector grid function on this fec and test setting/getting values
89 int vdim = mesh.SpaceDimension();
90 FiniteElementSpace vectorFes(&mesh, &fec, vdim);
91 GridFunction vectorGF(&vectorFes);
92 {
93 GridFunction& gf = vectorGF;
94 REQUIRE( gf.FESpace()->GetVDim() == vdim );
95 REQUIRE( gf.Size() == expScalarDofs * vdim );
96 gf = 0.;
97
98 // Test setting some vector values
99 mfem::Array<int> vdofs;
100 gf.FESpace()->GetElementVDofs(eltId,vdofs);
101 REQUIRE( vdofs.Size() > 0 );
102
103 // If we only expect dofs at one point, set non-zero as in the scalar case
104 if (expScalarDofs == 1)
105 {
106 for (int i = 0; i < vdofs.Size(); ++i)
107 {
108 gf[ vdofs[i] ] = (1. + i)/(1. + vdofs.Size());
109 }
110
111 // Check that the vector is not zero at the element center
112 Vector centerVec;
113 gf.GetVectorValue( eltId, Geometries.GetCenter(geom), centerVec);
114 REQUIRE( centerVec.Norml2() > 0. );
115 }
116 // Otherwise create a vector from each vertex to the element center
117 else
118 {
119 // Get position of element center
120 Vector ctr;
121 IsoparametricTransformation tr;
122 mesh.GetElementTransformation(eltId, &tr);
123 tr.Transform(Geometries.GetCenter(geom), ctr);
124
125 // Get each vertex position and subtract from ctr
126 Array<int> vertIndices;
127 mesh.GetElementVertices(eltId, vertIndices);
128 const int nv = vertIndices.Size();
129 for (int i=0; i < nv; ++i)
130 {
131 double* vertPos = mesh.GetVertex(vertIndices[i]);
132
133 for (int j=0; j < vdim; ++j)
134 {
135 gf[ vdofs[j*nv + i] ] = ctr[j] - vertPos[j];
136 }
137 }
138
139 // Check that none of the vectors are zero at vertices
140 Vector vec;
141 DenseMatrix dm;
142 gf.GetVectorValues(tr, *Geometries.GetVertices(geom), dm);
143 for (int i=0; i< nv; ++i)
144 {
145 dm.GetColumnReference(i, vec);
146 REQUIRE( vec.Norml2() > 0. );
147 }
148
149 // But, the vectors should cancel each other out at the midpoint
150 gf.GetVectorValue( eltId, Geometries.GetCenter(geom), vec);
151 REQUIRE( vec.Norml2() == MFEM_Approx(0.0) );
152 }
153
154 }
155
156 if (dumpMeshes)
157 {
158 const std::string prefix_path = "output_meshes/";
159
160 // Save mfem meshes using VisitDataCollection
161 VisItDataCollection dc(generateMeshFilename(fec, mesh), & mesh);
162 dc.SetPrefixPath(prefix_path);
163 dc.RegisterField("scalar_gf", &scalarGF);
164 dc.RegisterField("vector_gf", &vectorGF);
165 dc.Save();
166
167 // Save meshes and grid functions in VTK format
168 std::string fname = generateMeshFilename(fec, mesh, prefix_path, ".vtk");
169 std::fstream vtkFs( fname.c_str(), std::ios::out);
170
171 const int ref = 0;
172 mesh.PrintVTK( vtkFs, ref);
173 scalarGF.SaveVTK( vtkFs, "scalar_gf", ref);
174 vectorGF.SaveVTK( vtkFs, "vector_gf", ref);
175 }
176 }
177
178
179 TEST_CASE("Point mesh with linear grid function",
180 "[LinearFECollection]"
181 "[H1_FECollection]"
182 "[L2_FECollection]"
183 "[FiniteElementSpace]"
184 "[GridFunction]")
185 {
186 // Create a simple one element point mesh
187 int dim = 0, nv = 1, ne = 1, nb = 0, sdim = 2, order =1, attrib = 1;
188 Mesh mesh(dim, nv, ne, nb, sdim);
189
190 mesh.AddVertex(Vertex(0.,0.)());
191
192 int idx[1] = {0};
193 mesh.AddElement(new mfem::Point(idx, attrib));
194
195 mesh.FinalizeMesh();
196
197 REQUIRE( mesh.GetNV() == nv);
198 REQUIRE( mesh.GetNE() == ne);
199 REQUIRE( mesh.Dimension() == dim);
200 REQUIRE( mesh.SpaceDimension() == sdim);
201
202 SECTION("GridFunction with LinearFECollection")
203 {
204 LinearFECollection fec;
205 testGridFunctions(fec, mesh, nv);
206 }
207
208 SECTION("GridFunction with order 1 H1_FECollection")
209 {
210 H1_FECollection fec(order,dim);
211 testGridFunctions(fec, mesh, nv);
212 }
213
214 SECTION("GridFunction with L2_FECollection")
215 {
216 L2_FECollection fec(0,dim);
217 testGridFunctions(fec, mesh, 1);
218 }
219 }
220
221
222 TEST_CASE("Segment mesh with linear grid function",
223 "[LinearFECollection]"
224 "[H1_FECollection]"
225 "[L2_FECollection]"
226 "[FiniteElementSpace]"
227 "[GridFunction]")
228 {
229 // Create a simple one element segment mesh
230 int dim = 1, nv = 2, ne = 1, nb = 0, sdim = 2, order =1, attrib = 1;
231 Mesh mesh(dim, nv, ne, nb, sdim);
232
233 mesh.AddVertex(Vertex(0.,0.)());
234 mesh.AddVertex(Vertex(1.,0.)());
235
236 int idx[2] = {0,1};
237 mesh.AddElement(new Segment(idx, attrib));
238
239 mesh.FinalizeMesh();
240
241 REQUIRE( mesh.GetNV() == nv);
242 REQUIRE( mesh.GetNE() == ne);
243 REQUIRE( mesh.Dimension() == dim);
244 REQUIRE( mesh.SpaceDimension() == sdim);
245
246 SECTION("GridFunction with LinearFECollection")
247 {
248 LinearFECollection fec;
249 testGridFunctions(fec, mesh, nv);
250 }
251
252 SECTION("GridFunction with order 1 H1_FECollection")
253 {
254 H1_FECollection fec(order,dim);
255 testGridFunctions(fec, mesh, nv);
256 }
257
258 SECTION("GridFunction with L2_FECollection")
259 {
260 L2_FECollection fec(0,dim);
261 testGridFunctions(fec, mesh, 1);
262 }
263 }
264
265 TEST_CASE("Triangle mesh with linear grid function",
266 "[LinearFECollection]"
267 "[H1_FECollection]"
268 "[L2_FECollection]"
269 "[FiniteElementSpace]"
270 "[GridFunction]")
271 {
272 // Create a simple one element triangle mesh
273 int dim = 2, nv = 3, ne = 1, nb = 0, sdim = 2, order =1, attrib = 1;
274 Mesh mesh(dim, nv, ne, nb, sdim);
275
276 mesh.AddVertex(Vertex(0.,0.)());
277 mesh.AddVertex(Vertex(1.,0.)());
278 mesh.AddVertex(Vertex(1.,1.)());
279
280 int idx[3] = {0,1,2};
281 mesh.AddElement(new Triangle(idx, attrib));
282
283 mesh.FinalizeMesh();
284
285 REQUIRE( mesh.GetNV() == nv);
286 REQUIRE( mesh.GetNE() == ne);
287 REQUIRE( mesh.Dimension() == dim);
288 REQUIRE( mesh.SpaceDimension() == sdim);
289
290 SECTION("GridFunction with LinearFECollection")
291 {
292 LinearFECollection fec;
293 testGridFunctions(fec, mesh, nv);
294 }
295
296 SECTION("GridFunction with order 1 H1_FECollection")
297 {
298 H1_FECollection fec(order,dim);
299 testGridFunctions(fec, mesh, nv);
300 }
301
302 SECTION("GridFunction with L2_FECollection")
303 {
304 L2_FECollection fec(0,dim);
305 testGridFunctions(fec, mesh, 1);
306 }
307 }
308
309 TEST_CASE("Quad mesh with linear grid function",
310 "[LinearFECollection]"
311 "[H1_FECollection]"
312 "[L2_FECollection]"
313 "[FiniteElementSpace]"
314 "[GridFunction]")
315 {
316 // Create a simple one element quad mesh
317 int dim = 2, nv = 4, ne = 1, nb = 0, sdim = 2, order =1, attrib = 1;
318 Mesh mesh(dim, nv, ne, nb, sdim);
319
320 mesh.AddVertex(Vertex(0.,0.)());
321 mesh.AddVertex(Vertex(1.,0.)());
322 mesh.AddVertex(Vertex(1.,1.)());
323 mesh.AddVertex(Vertex(0.,1.)());
324
325 int idx[4] = {0,1,2,3};
326 mesh.AddElement(new Quadrilateral(idx, attrib));
327
328 mesh.FinalizeMesh();
329
330 REQUIRE( mesh.GetNV() == nv);
331 REQUIRE( mesh.GetNE() == ne);
332 REQUIRE( mesh.Dimension() == dim);
333 REQUIRE( mesh.SpaceDimension() == sdim);
334
335 SECTION("GridFunction with LinearFECollection")
336 {
337 LinearFECollection fec;
338 testGridFunctions(fec, mesh, nv);
339 }
340
341 SECTION("GridFunction with order 1 H1_FECollection")
342 {
343 H1_FECollection fec(order,dim);
344 testGridFunctions(fec, mesh, nv);
345 }
346
347 SECTION("GridFunction with L2_FECollection")
348 {
349 L2_FECollection fec(0,dim);
350 testGridFunctions(fec, mesh, 1);
351 }
352
353 }
354
355
356 TEST_CASE("Tet mesh with linear grid function",
357 "[LinearFECollection]"
358 "[H1_FECollection]"
359 "[L2_FECollection]"
360 "[FiniteElementSpace]"
361 "[GridFunction]")
362 {
363 // Create a simple one element tet mesh
364 int dim = 3, nv = 4, ne = 1, nb = 0, sdim = 3, order =1, attrib = 1;
365 Mesh mesh(dim, nv, ne, nb, sdim);
366
367 mesh.AddVertex(Vertex(0.,0.,0.)());
368 mesh.AddVertex(Vertex(1.,0.,0.)());
369 mesh.AddVertex(Vertex(1.,1.,0.)());
370 mesh.AddVertex(Vertex(1.,1.,1.)());
371
372 int idx[4] = {0,1,2,3};
373 mesh.AddElement(new Tetrahedron(idx, attrib));
374
375 mesh.FinalizeMesh();
376
377 REQUIRE( mesh.GetNV() == nv);
378 REQUIRE( mesh.GetNE() == ne);
379 REQUIRE( mesh.Dimension() == dim);
380 REQUIRE( mesh.SpaceDimension() == sdim);
381
382 SECTION("GridFunction with LinearFECollection")
383 {
384 LinearFECollection fec;
385 testGridFunctions(fec, mesh, nv);
386 }
387
388 SECTION("GridFunction with order 1 H1_FECollection")
389 {
390 H1_FECollection fec(order,dim);
391 testGridFunctions(fec, mesh, nv);
392 }
393
394 SECTION("GridFunction with L2_FECollection")
395 {
396 L2_FECollection fec(0,dim);
397 testGridFunctions(fec, mesh, 1);
398 }
399
400 }
401
402 TEST_CASE("Hex mesh with linear grid function",
403 "[LinearFECollection]"
404 "[H1_FECollection]"
405 "[L2_FECollection]"
406 "[FiniteElementSpace]"
407 "[GridFunction]")
408 {
409 // Create a simple one element hex mesh
410 int dim = 3, nv = 8, ne = 1, nb = 0, sdim = 3, order =1, attrib = 1;
411 Mesh mesh(dim, nv, ne, nb, sdim);
412
413 mesh.AddVertex(Vertex(0.,0.,0.)());
414 mesh.AddVertex(Vertex(1.,0.,0.)());
415 mesh.AddVertex(Vertex(1.,1.,0.)());
416 mesh.AddVertex(Vertex(0.,1.,0.)());
417 mesh.AddVertex(Vertex(0.,0.,1.)());
418 mesh.AddVertex(Vertex(1.,0.,1.)());
419 mesh.AddVertex(Vertex(1.,1.,1.)());
420 mesh.AddVertex(Vertex(0.,1.,1.)());
421
422 int idx[8] = {0,1,2,3,4,5,6,7};
423 mesh.AddElement(new Hexahedron(idx, attrib));
424
425 mesh.FinalizeMesh();
426
427 REQUIRE( mesh.GetNV() == nv);
428 REQUIRE( mesh.GetNE() == ne);
429 REQUIRE( mesh.Dimension() == dim);
430 REQUIRE( mesh.SpaceDimension() == sdim);
431
432 SECTION("GridFunction with LinearFECollection")
433 {
434 LinearFECollection fec;
435 testGridFunctions(fec, mesh, nv);
436 }
437
438 SECTION("GridFunction with order 1 H1_FECollection")
439 {
440 H1_FECollection fec(order,dim);
441 testGridFunctions(fec, mesh, nv);
442 }
443
444 SECTION("GridFunction with L2_FECollection")
445 {
446 L2_FECollection fec(0,dim);
447 testGridFunctions(fec, mesh, 1);
448 }
449
450 }
451