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 
15 using namespace mfem;
16 
17 namespace qf_coeff
18 {
19 
20 TEST_CASE("Quadrature Function Coefficients",
21           "[Quadrature Function Coefficients]")
22 {
23    int order_h1 = 2, n = 4, dim = 3;
24    double tol = 1e-14;
25 
26    Mesh mesh = Mesh::MakeCartesian3D(
27                   n, n, n, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
28    mesh.SetCurvature(order_h1);
29 
30    int intOrder = 2 * order_h1 + 1;
31 
32    QuadratureSpace qspace(&mesh, intOrder);
33    QuadratureFunction quadf_coeff(&qspace, 1);
34    QuadratureFunction quadf_vcoeff(&qspace, dim);
35 
36    const IntegrationRule ir = qspace.GetElementIntRule(0);
37 
38    const GeometricFactors *geom_facts =
39       mesh.GetGeometricFactors(ir, GeometricFactors::COORDINATES);
40 
41    {
42       int nelems = quadf_coeff.Size() / quadf_coeff.GetVDim() / ir.GetNPoints();
43       int vdim = ir.GetNPoints();
44 
45       for (int i = 0; i < nelems; i++)
46       {
47          for (int j = 0; j < vdim; j++)
48          {
49             //X has dims nqpts x sdim x ne
50             quadf_coeff((i * vdim) + j) =
51                geom_facts->X((i * vdim * dim) + (vdim * 2) + j );
52          }
53       }
54    }
55 
56    {
57       int nqpts = ir.GetNPoints();
58       int nelems = quadf_vcoeff.Size() / quadf_vcoeff.GetVDim() / nqpts;
59       int vdim = quadf_vcoeff.GetVDim();
60 
61       for (int i = 0; i < nelems; i++)
62       {
63          for (int j = 0; j < vdim; j++)
64          {
65             for (int k = 0; k < nqpts; k++)
66             {
67                //X has dims nqpts x sdim x ne
68                quadf_vcoeff((i * nqpts * vdim) + (k * vdim ) + j) =
69                   geom_facts->X((i * nqpts * vdim) + (j * nqpts) + k);
70             }
71          }
72       }
73    }
74 
75    QuadratureFunctionCoefficient qfc(quadf_coeff);
76    VectorQuadratureFunctionCoefficient qfvc(quadf_vcoeff);
77 
78    SECTION("Operators on VecQuadFuncCoeff")
79    {
80       std::cout << "Testing VecQuadFuncCoeff: " << std::endl;
81 #ifdef MFEM_USE_EXCEPTIONS
82       std::cout << " Setting Component" << std::endl;
83       REQUIRE_THROWS(qfvc.SetComponent(3, 1));
84       REQUIRE_THROWS(qfvc.SetComponent(-1, 1));
85       REQUIRE_NOTHROW(qfvc.SetComponent(1, 2));
86       REQUIRE_THROWS(qfvc.SetComponent(0, 4));
87       REQUIRE_THROWS(qfvc.SetComponent(1, 3));
88       REQUIRE_NOTHROW(qfvc.SetComponent(0, 2));
89       REQUIRE_THROWS(qfvc.SetComponent(0, 0));
90 #endif
91       qfvc.SetComponent(0, 3);
92    }
93 
94    SECTION("Operators on VectorQuadratureLFIntegrator")
95    {
96       std::cout << "Testing VectorQuadratureLFIntegrator: " << std::endl;
97       H1_FECollection    fec_h1(order_h1, dim);
98       FiniteElementSpace fespace_h1(&mesh, &fec_h1, dim);
99 
100       GridFunction nodes(&fespace_h1);
101       mesh.GetNodes(nodes);
102 
103       Vector output(nodes.Size());
104       output = 0.0;
105 
106       LinearForm lf(&fespace_h1);
107       lf.AddDomainIntegrator(new VectorQuadratureLFIntegrator(qfvc, NULL));
108 
109       lf.Assemble();
110 
111       BilinearForm L2(&fespace_h1);
112 
113       L2.AddDomainIntegrator(new VectorMassIntegrator());
114       L2.Assemble();
115 
116       SparseMatrix mat = L2.SpMat();
117 
118       mat.Mult(nodes, output);
119 
120       output -= lf;
121 
122       REQUIRE(output.Norml2() < tol);
123    }
124 
125    SECTION("Operators on QuadratureLFIntegrator")
126    {
127       std::cout << "Testing QuadratureLFIntegrator: " << std::endl;
128       H1_FECollection    fec_h1(order_h1, dim);
129       FiniteElementSpace fespace_h1(&mesh, &fec_h1, 1);
130       FiniteElementSpace fespace_h3(&mesh, &fec_h1, 3);
131 
132       GridFunction nodes(&fespace_h3);
133       mesh.GetNodes(nodes);
134 
135       Vector output(nodes.Size() / dim);
136       Vector nz(nodes.Size() / dim);
137       output = 0.0;
138 
139       nz.MakeRef(nodes, nz.Size() * 2);
140 
141       LinearForm lf(&fespace_h1);
142       lf.AddDomainIntegrator(new QuadratureLFIntegrator(qfc, NULL));
143 
144       lf.Assemble();
145 
146       BilinearForm L2(&fespace_h1);
147 
148       L2.AddDomainIntegrator(new MassIntegrator(&ir));
149       L2.Assemble();
150 
151       SparseMatrix mat = L2.SpMat();
152 
153       mat.Mult(nz, output);
154 
155       output -= lf;
156 
157       REQUIRE(output.Norml2() < tol);
158    }
159 
160 }
161 
162 } // namespace qf_coeff
163