1 /** @file gsPoissonSolver_test.cpp
2
3 @brief Tests for Poisson solvers.
4
5 Verification test for the vector valued Poisson solver
6 with the gsPdeAssembler structure.
7 The test uses the Method of Manufactured Solutions (MMS)
8 and measure the convergence rate of the L2 error. If the
9 convergence rate is not sufficient, the test fails.
10
11 Test inclueds:
12
13 Inhomogeneous source term (right hand side)
14 Inhomogeneous Dirichlet boundary conditions
15 Inhomogeneous Neumann boundary conditions
16 Nitche handling of Dirichlet BC
17 DG handling of multipatch
18
19 x,y \in (0,1)
20
21 Source function:
22
23 f(x,y)_x = ((pi*k0)^2 + (pi*k1)^2)*sin(pi*x*k0)*sin(pi*y*k1)
24 f(x,y)_y = ((pi*k2)^2 + (pi*k3)^2)*sin(pi*x*k2)*sin(pi*y*k3)
25
26 Solution:
27
28 u(x,y)_x = sin(pi*x*k0)*sin(pi*y*k1)+pi/10
29 u(x,y)_y = sin(pi*x*k2)*sin(pi*y*k3)-pi/10
30
31 This Source Code Form is subject to the terms of the Mozilla Public
32 License, v. 2.0. If a copy of the MPL was not distributed with this
33 file, You can obtain one at http://mozilla.org/MPL/2.0/.
34
35 Author(s): J. Sogn, S. Takacs
36 */
37
38 #include "gismo_unittest.h"
39 #include <gsSolver/gsSolverUtils.h>
40
41 using namespace gismo;
42
43
runPoissonSolverTest(dirichlet::strategy Dstrategy,iFace::strategy Istrategy)44 void runPoissonSolverTest( dirichlet::strategy Dstrategy, iFace::strategy Istrategy )
45 {
46 int numRefine = 2;
47 int maxIterations = 3;
48 real_t convratetol = 1.85; // Convergence rate should be around 2
49
50 // List of element sizes
51 std::vector <real_t> h_list;
52
53 // List of L2-errors
54 std::vector <real_t> error_list;
55
56 // Source function
57 gsFunctionExpr<> f("((pi*1)^2 + (pi*2)^2)*sin(pi*x*1)*sin(pi*y*2)",
58 "((pi*3)^2 + (pi*4)^2)*sin(pi*x*3)*sin(pi*y*4)",2);
59
60 // Exact solution
61 gsFunctionExpr<> g("sin(pi*x*1)*sin(pi*y*2)+pi/10",
62 "sin(pi*x*3)*sin(pi*y*4)-pi/10",2);
63
64 // Define Geometry (Unit square with 4 patches)
65 gsMultiPatch<> patches = gsNurbsCreator<>::BSplineSquareGrid(2, 2, 0.5);
66
67 // Define Boundary conditions
68 gsBoundaryConditions<> bcInfo;
69
70 //Dirichlet BCs
71 bcInfo.addCondition(0, boundary::west, condition_type::dirichlet, &g);
72 bcInfo.addCondition(1, boundary::west, condition_type::dirichlet, &g);
73 bcInfo.addCondition(1, boundary::north, condition_type::dirichlet, &g);
74 bcInfo.addCondition(3, boundary::north, condition_type::dirichlet, &g);
75
76 // Neumann BCs
77 gsFunctionExpr<> gEast ("1*pi*cos(pi*1)*sin(pi*2*y)", "3*pi*cos(pi*3)*sin(pi*4*y)",2);
78 gsFunctionExpr<> gSouth("-pi*2*sin(pi*x*1)","-pi*4*sin(pi*x*3)",2);
79
80 bcInfo.addCondition(3, boundary::east, condition_type::neumann, &gEast);
81 bcInfo.addCondition(2, boundary::east, condition_type::neumann, &gEast);
82 bcInfo.addCondition(0, boundary::south, condition_type::neumann, &gSouth);
83 bcInfo.addCondition(2, boundary::south, condition_type::neumann, &gSouth);
84
85 // Copy bases for refinement
86 gsMultiBasis<> refine_bases( patches );
87
88 // Define discretization space by initial refining the basis of the geometry
89 for (int i = 0; i < numRefine; ++i)
90 refine_bases.uniformRefine();
91
92
93 // linear solver
94 gsSparseSolver<>::CGDiagonal solver;
95 gsMatrix<> solVector;
96
97 // Start Loop
98 // For each iteration we h-refine the basis, then set up and solve the
99 // Poisson problem. Then find the L2 error and element size.
100 for (int i = 0; i < maxIterations; ++i)
101 {
102 refine_bases.uniformRefine();
103
104 // Initilize Assembler
105 gsPoissonAssembler<real_t> poisson(patches,refine_bases,bcInfo,f,Dstrategy,Istrategy);
106 //gsPoissonAssembler<> poisson(*patches, bcInfo, refine_bases, f);
107
108 // Assemble and solve
109 poisson.assemble();
110 solVector = solver.compute( poisson.matrix() ).solve( poisson.rhs() );
111
112 // Find the element size
113 real_t h = math::pow( (real_t) refine_bases.size(0), -1.0 / refine_bases.dim() );
114 h_list.push_back(h);
115
116 // Access the solution
117 const gsField<> sol = poisson.constructSolution(solVector);
118
119 // Find the l2 error
120 real_t l2error = sol.distanceL2(g);
121 error_list.push_back(l2error);
122
123 }
124
125 // Finding convergence rate in to differnt ways
126
127 // Convergence rate found by last to errors are elemet sizes
128 real_t convratelast = math::log(error_list[error_list.size()-2]/error_list[error_list.size()-1])/
129 math::log(h_list[h_list.size()-2]/h_list[h_list.size()-1]);
130
131 // Convergence rate found by a least square method
132 // PS: This method of finding convergence rate mights require some initialt
133 // refinement such that convergence have startet for the largest element
134 // size value.
135 real_t convrateavg = gsSolverUtils<>::convergenceRateLS(error_list,h_list);
136
137 CHECK( convrateavg > convratetol );
138 CHECK( convratelast > convratetol );
139
140 //Clean up before next test
141 error_list.clear();
142 h_list.clear();
143 }
144
145
SUITE(gsPoissonSolver_test)146 SUITE(gsPoissonSolver_test)
147 {
148
149 TEST(Galerkin_test)
150 {
151 runPoissonSolverTest(dirichlet::elimination, iFace::glue);
152 }
153
154 TEST(dG_test)
155 {
156 runPoissonSolverTest(dirichlet::elimination, iFace::dg);
157 }
158
159 TEST(Nitsche_dG_test)
160 {
161 runPoissonSolverTest(dirichlet::nitsche, iFace::dg);
162 }
163
164 }
165
166