1 /** @file gsElPoissonAssembler.hpp
2
3 @brief Provides stiffness matrix for Poisson's equation.
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s):
12 A.Shamanskiy (2016 - ...., TU Kaiserslautern)
13 */
14
15 #pragma once
16
17 #include <gsElasticity/gsElPoissonAssembler.h>
18
19 #include <gsPde/gsPoissonPde.h>
20 #include <gsElasticity/gsVisitorElPoisson.h>
21
22 namespace gismo
23 {
24
25 template<class T>
gsElPoissonAssembler(const gsMultiPatch<T> & patches,const gsMultiBasis<T> & basis,const gsBoundaryConditions<T> & bconditions,const gsFunction<T> & body_force)26 gsElPoissonAssembler<T>::gsElPoissonAssembler(const gsMultiPatch<T> & patches,
27 const gsMultiBasis<T> & basis,
28 const gsBoundaryConditions<T> & bconditions,
29 const gsFunction<T> & body_force)
30 {
31 gsPiecewiseFunction<T> rightHandSides;
32 rightHandSides.addPiece(body_force);
33 typename gsPde<T>::Ptr pde( new gsPoissonPde<T>(patches,bconditions,rightHandSides) );
34 m_bases.push_back(basis);
35 Base::initialize(pde, m_bases, defaultOptions());
36 }
37
38 template <class T>
defaultOptions()39 gsOptionList gsElPoissonAssembler<T>::defaultOptions()
40 {
41 gsOptionList opt = Base::defaultOptions();
42 opt.addReal("LocalStiff","Stiffening degree for the Jacobian-based local stiffening",0.);
43 return opt;
44 }
45
46 template <class T>
refresh()47 void gsElPoissonAssembler<T>::refresh()
48 {
49 std::vector<gsDofMapper> m_dofMappers(m_bases.size());
50 m_bases[0].getMapper((dirichlet::strategy)m_options.getInt("DirichletStrategy"),
51 iFace::glue,m_pde_ptr->bc(),m_dofMappers[0],0,true);
52
53 m_system = gsSparseSystem<T>(m_dofMappers[0]);
54 m_system.reserve(m_bases[0], m_options, m_pde_ptr->numRhs());
55 Base::computeDirichletDofs(0);
56 }
57
58 template<class T>
assemble(bool saveEliminationMatrix)59 void gsElPoissonAssembler<T>::assemble(bool saveEliminationMatrix)
60 {
61 m_system.matrix().setZero();
62 m_system.reserve(m_bases[0], m_options, m_pde_ptr->numRhs());
63 m_system.rhs().setZero(Base::numDofs(),m_pde_ptr->numRhs());
64
65 if (saveEliminationMatrix)
66 {
67 eliminationMatrix.resize(Base::numDofs(),Base::numFixedDofs());
68 eliminationMatrix.setZero();
69 eliminationMatrix.reservePerColumn(m_system.numColNz(m_bases[0],m_options));
70 }
71
72 gsVisitorElPoisson<T> visitor(*m_pde_ptr, saveEliminationMatrix ? &eliminationMatrix : nullptr);
73 Base::template push<gsVisitorElPoisson<T> >(visitor);
74
75 m_system.matrix().makeCompressed();
76
77 if (saveEliminationMatrix)
78 {
79 Base::rhsWithZeroDDofs = m_system.rhs();
80 eliminationMatrix.makeCompressed();
81 }
82 }
83
84 template <class T>
constructSolution(const gsMatrix<T> & solVector,const std::vector<gsMatrix<T>> & fixedDoFs,gsMultiPatch<T> & result) const85 void gsElPoissonAssembler<T>::constructSolution(const gsMatrix<T> & solVector,
86 const std::vector<gsMatrix<T> > & fixedDoFs,
87 gsMultiPatch<T> & result) const
88 {
89 Base::constructSolution(solVector,fixedDoFs,result,gsVector<index_t>::Zero(1));
90 }
91
92 }// namespace gismo ends
93