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 TEST_CASE("OperatorChebyshevSmoother", "[Chebyshev symmetry]") 18 { 19 for (int order = 2; order < 5; ++order) 20 { 21 const int cheb_order = 2; 22 23 Mesh mesh = Mesh::MakeCartesian3D(4, 4, 4, Element::HEXAHEDRON); 24 FiniteElementCollection *fec = new H1_FECollection(order, 3); 25 FiniteElementSpace fespace(&mesh, fec); 26 Array<int> ess_bdr(mesh.bdr_attributes.Max()); 27 ess_bdr = 1; 28 Array<int> ess_tdof_list; 29 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list); 30 31 BilinearForm aform(&fespace); 32 aform.SetAssemblyLevel(AssemblyLevel::PARTIAL); 33 aform.AddDomainIntegrator(new DiffusionIntegrator); 34 aform.Assemble(); 35 OperatorPtr opr; 36 opr.SetType(Operator::ANY_TYPE); 37 aform.FormSystemMatrix(ess_tdof_list, opr); 38 Vector diag(fespace.GetTrueVSize()); 39 aform.AssembleDiagonal(diag); 40 41 Solver* smoother = new OperatorChebyshevSmoother(*opr, diag, ess_tdof_list, 42 cheb_order); 43 44 int n = smoother->Width(); 45 Vector left(n); 46 Vector right(n); 47 int seed = (int) time(0); 48 left.Randomize(seed); 49 right.Randomize(seed + 2); 50 51 // test that x^T S y = y^T S x 52 Vector out(n); 53 out = 0.0; 54 smoother->Mult(right, out); 55 double forward_val = left * out; 56 smoother->Mult(left, out); 57 double transpose_val = right * out; 58 59 double error = fabs(forward_val - transpose_val) / fabs(forward_val); 60 std::cout << "Order " << order << " symmetry error: " << error << std::endl; 61 REQUIRE(error < 1.e-13); 62 63 delete smoother; 64 } 65 } 66