/
gsBiharmonicAssembler.hpp
131 lines (107 loc) · 4.47 KB
/
gsBiharmonicAssembler.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
/** @file gsElMassAssembler.hpp
@brief Provides stiffness matrix for Poisson's equation.
This file is part of the G+Smo library.
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
Author(s):
A.Shamanskiy (2016 - ...., TU Kaiserslautern)
*/
#pragma once
#include <gsElasticity/gsBiharmonicAssembler.h>
#include <gsPde/gsPoissonPde.h>
#include <gsElasticity/gsVisitorBiharmonic.h>
namespace gismo
{
template<class T>
gsBiharmonicAssembler<T>::gsBiharmonicAssembler(const gsMultiPatch<T> & patches,
const gsMultiBasis<T> & basis,
const gsBoundaryConditions<T> & bconditions,
const gsFunction<T> & body_force)
{
gsPiecewiseFunction<T> rightHandSides;
rightHandSides.addPiece(body_force);
typename gsPde<T>::Ptr pde( new gsPoissonPde<T>(patches,bconditions,rightHandSides) );
m_bases.push_back(basis);
m_bases.push_back(basis);
Base::initialize(pde, m_bases, defaultOptions());
}
template <class T>
gsOptionList gsBiharmonicAssembler<T>::defaultOptions()
{
gsOptionList opt = Base::defaultOptions();
opt.addReal("LocalStiff","Stiffening degree for the Jacobian-based local stiffening",0.);
return opt;
}
template <class T>
void gsBiharmonicAssembler<T>::refresh()
{
std::vector<gsDofMapper> m_dofMappers(2);
for (unsigned d = 0; d < 2; d++)
m_bases[d].getMapper((dirichlet::strategy)m_options.getInt("DirichletStrategy"),
iFace::glue,m_pde_ptr->bc(),m_dofMappers[d],d,true);
m_system = gsSparseSystem<T>(m_dofMappers, gsVector<index_t>::Ones(2));
reserve();
Base::computeDirichletDofs(0);
Base::computeDirichletDofs(1);
}
template <class T>
void gsBiharmonicAssembler<T>::reserve()
{
// Pick up values from options
const T bdA = m_options.getReal("bdA");
const index_t bdB = m_options.getInt("bdB");
const T bdO = m_options.getReal("bdO");
index_t dim = m_bases[0][0].dim();
index_t deg = 0;
for (index_t d = 0; d < dim; ++d )
if (m_bases[0][0].degree(d) > deg)
deg = m_bases[0][0].degree(d);
index_t numElPerColumn = pow((bdA*deg+bdB),dim)*2;
m_system.reserve(numElPerColumn*(1+bdO),m_pde_ptr->numRhs());
}
template<class T>
void gsBiharmonicAssembler<T>::assemble(bool saveEliminationMatrix)
{
m_system.matrix().setZero();
reserve();
m_system.rhs().setZero(Base::numDofs(),m_pde_ptr->numRhs());
if (saveEliminationMatrix)
{
eliminationMatrix.resize(Base::numDofs(),Base::numFixedDofs());
eliminationMatrix.setZero();
eliminationMatrix.reservePerColumn(m_system.numColNz(m_bases[0],m_options));
}
gsVisitorBiharmonic<T> visitor(*m_pde_ptr, saveEliminationMatrix ? &eliminationMatrix : nullptr);
Base::template push<gsVisitorBiharmonic<T> >(visitor);
m_system.matrix().makeCompressed();
if (saveEliminationMatrix)
{
Base::rhsWithZeroDDofs = m_system.rhs();
eliminationMatrix.makeCompressed();
}
}
//--------------------- SOLUTION CONSTRUCTION ----------------------------------//
template <class T>
void gsBiharmonicAssembler<T>::constructSolution(const gsMatrix<T> & solVector,
const std::vector<gsMatrix<T> > & fixedDoFs,
gsMultiPatch<T> & solution) const
{
Base::constructSolution(solVector,fixedDoFs,solution,gsVector<index_t>::Zero(1));
}
template <class T>
void gsBiharmonicAssembler<T>::constructSolutionAux(const gsMatrix<T> & solVector,
const std::vector<gsMatrix<T> > & fixedDoFs,
gsMultiPatch<T> & solutionAux) const
{
Base::constructSolution(solVector,fixedDoFs,solutionAux,gsVector<index_t>::Ones(1));
}
template <class T>
void gsBiharmonicAssembler<T>::constructSolution(const gsMatrix<T> & solVector,
const std::vector<gsMatrix<T> > & fixedDoFs,
gsMultiPatch<T> & solutionMain, gsMultiPatch<T> & solutionAux) const
{
constructSolution(solVector,fixedDoFs,solutionMain);
constructSolutionAux(solVector,fixedDoFs,solutionAux);
}
}// namespace gismo ends