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