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