1 /** @file poisson_example.cpp
2
3 @brief Tutorial on how to use G+Smo to solve the Poisson equation,
4 see the \ref PoissonTutorial
5
6 This file is part of the G+Smo library.
7
8 This Source Code Form is subject to the terms of the Mozilla Public
9 License, v. 2.0. If a copy of the MPL was not distributed with this
10 file, You can obtain one at http://mozilla.org/MPL/2.0/.
11
12 Author(s): J. Sogn
13 */
14
15 //! [Include namespace]
16 # include <gismo.h>
17
18 using namespace gismo;
19 //! [Include namespace]
20
main(int argc,char * argv[])21 int main(int argc, char *argv[])
22 {
23 //! [Parse command line]
24 bool plot = false;
25
26 gsCmdLine cmd("Tutorial on solving a Poisson problem.");
27 cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot);
28 try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
29 //! [Parse command line]
30
31 //! [Function data]
32 // Define source function
33 gsFunctionExpr<> f("((pi*1)^2 + (pi*2)^2)*sin(pi*x*1)*sin(pi*y*2)",
34 "((pi*3)^2 + (pi*4)^2)*sin(pi*x*3)*sin(pi*y*4)",2);
35 // For homogeneous term, we can use this (last argument is the dimension of the domain)
36 //gsConstantFunction<> f(0.0, 0.0, 2);
37
38 // Define exact solution (optional)
39 gsFunctionExpr<> g("sin(pi*x*1)*sin(pi*y*2)+pi/10",
40 "sin(pi*x*3)*sin(pi*y*4)-pi/10",2);
41
42 // Print out source function and solution
43 gsInfo<<"Source function "<< f << "\n";
44 gsInfo<<"Exact solution "<< g <<"\n\n";
45 //! [Function data]
46
47 //! [Geometry data]
48 // Define Geometry, must be a gsMultiPatch object
49 gsMultiPatch<> patches;
50 // Create 4 (2 x 2) patches of squares:
51 //
52 // Square/patch 0 is in lower left corner
53 // Square/patch 1 is in upper left corner
54 // Square/patch 2 is in lower right corner
55 // Square/patch 3 is in upper right corner
56 //
57 // The last argument scale the squares such that we
58 // get the unit square as domain.
59 patches = gsNurbsCreator<>::BSplineSquareGrid(2, 2, 0.5);
60 gsInfo << "The domain is a "<< patches <<"\n";
61 //! [Geometry data]
62
63 // For single patch unit square of quadratic elements use (Note:
64 // you need to update the bounadry conditions section for this to
65 // work properly!) :
66 // patches = gsMultiPatch<>(*gsNurbsCreator<>::BSplineSquare(2));
67
68 // Geometry can also be read from file (if gsMultiPatch):
69 // gsReadFile<>("planar/lshape_p2.xml", patches);
70
71 // Define Boundary conditions. Note that if one boundary is
72 // "free", eg. if no condition is defined, then it is a natural
73 // boundary (zero Neumann condition)
74 // Also, remember that a pure Neumann problem has no unique
75 // solution, thereforer implies a singular matrix. In this case
76 // a corner DoF can be fixed to a given value to obtain a unique solution.
77 // (example: bcInfo.addCornerValue(boundary::southwest, value, patch);)
78
79 //! [Boundary conditions]
80 gsBoundaryConditions<> bcInfo;
81 // Every patch with a boundary need to be specified. In this
82 // there are in total 8 sides (two for each patch)
83
84 // Dirichlet Boundary conditions
85 // First argument is the patch number
86 bcInfo.addCondition(0, boundary::west, condition_type::dirichlet, &g);
87 bcInfo.addCondition(1, boundary::west, condition_type::dirichlet, &g);
88
89 bcInfo.addCondition(1, boundary::north, condition_type::dirichlet, &g);
90 bcInfo.addCondition(3, boundary::north, condition_type::dirichlet, &g);
91
92 // Neumann Boundary conditions
93 gsFunctionExpr<> hEast ("1*pi*cos(pi*1)*sin(pi*2*y)", "3*pi*cos(pi*3)*sin(pi*4*y)",2);
94 gsFunctionExpr<> hSouth("-pi*2*sin(pi*x*1)","-pi*4*sin(pi*x*3)",2);
95
96 bcInfo.addCondition(3, boundary::east, condition_type::neumann, &hEast);
97 bcInfo.addCondition(2, boundary::east, condition_type::neumann, &hEast);
98
99 bcInfo.addCondition(0, boundary::south, condition_type::neumann, &hSouth);
100 bcInfo.addCondition(2, boundary::south, condition_type::neumann, &hSouth);
101 //! [Boundary conditions]
102
103 /*
104 //Alternatively: You can automatically create Dirichlet boundary
105 //conditions using one function (the exact solution) for all
106 //boundaries like this:
107
108 for (gsMultiPatch<>::const_biterator
109 bit = patches.bBegin(); bit != patches.bEnd(); ++bit)
110 {
111 bcInfo.addCondition( *bit, condition_type::dirichlet, &g );
112 }
113 */
114
115 //! [Refinement]
116 // Copy basis from the geometry
117 gsMultiBasis<> refine_bases( patches );
118
119 // Number for h-refinement of the computational (trial/test) basis.
120 const int numRefine = 2;
121
122 // Number for p-refinement of the computational (trial/test) basis.
123 const int degree = 2;
124
125 // h-refine each basis (4, one for each patch)
126 for ( int i = 0; i < numRefine; ++i)
127 refine_bases.uniformRefine();
128
129 // k-refinement (set degree)
130 for ( size_t i = 0; i < refine_bases.nBases(); ++ i )
131 refine_bases[i].setDegreePreservingMultiplicity(degree);
132
133 //! [Refinement]
134
135 ////////////// Setup solver and solve //////////////
136 // Initialize Solver
137 // Setup method for handling Dirichlet boundaries, options:
138 //
139 // * elimination: Eliminate the Dirichlet DoFs from the linear system.
140 //
141 // * nitsche: Keep the Dirichlet DoFs and enforce the boundary
142 //
143 // condition weakly by a penalty term.
144 // Setup method for handling patch interfaces, options:
145 //
146 // * glue:Glue patches together by merging DoFs across an interface into one.
147 // This only works for conforming interfaces
148 //
149 // * dg: Use discontinuous Galerkin-like coupling between adjacent patches.
150 // (This option might not be available yet)
151 //! [Assemble]
152 gsPoissonAssembler<real_t> assembler(patches,refine_bases,bcInfo,f,
153 //dirichlet::elimination, iFace::glue);
154 dirichlet::nitsche , iFace::glue);
155
156 // Generate system matrix and load vector
157 gsInfo<< "Assembling...\n";
158 assembler.assemble();
159 gsInfo << "Have assembled a system (matrix and load vector) with "
160 << assembler.numDofs() << " dofs.\n";
161 //! [Assemble]
162
163 //! [Solve]
164 // Initialize the conjugate gradient solver
165 gsInfo << "Solving...\n";
166 gsSparseSolver<>::CGDiagonal solver( assembler.matrix() );
167 gsMatrix<> solVector = solver.solve( assembler.rhs() );
168 gsInfo << "Solved the system with CG solver.\n";
169 //! [Solve]
170
171 //! [Construct solution]
172 // Construct the solution as a scalar field
173 gsMultiPatch<> mpsol;
174 assembler.constructSolution(solVector, mpsol);
175 gsField<> sol( assembler.patches(), mpsol);
176 //! [Construct solution]
177
178 if (plot)
179 {
180 //! [Plot in Paraview]
181 // Write approximate and exact solution to paraview files
182 gsInfo<<"Plotting in Paraview...\n";
183 gsWriteParaview<>(sol, "poisson2d", 1000);
184 const gsField<> exact( assembler.patches(), g, false );
185 gsWriteParaview<>( exact, "poisson2d_exact", 1000);
186
187 // Run paraview
188 gsFileManager::open("poisson2d.pvd");
189 gsFileManager::open("poisson2d_exact.pvd");
190 //! [Plot in Paraview]
191 }
192 else
193 {
194 gsInfo << "Done. No output created, re-run with --plot to get a ParaView "
195 "file containing the solution.\n";
196 }
197 return EXIT_SUCCESS;
198
199 }// end main
200