G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsL2Projection.hpp
Go to the documentation of this file.
1 
14 #include <gsPde/gsBoundaryConditions.h>
15 #include <gsAssembler/gsExprAssembler.h>
16 #include <gsAssembler/gsExprEvaluator.h>
17 #include <gsMatrix/gsSparseSolver.h>
18 #include <gsCore/gsBoundary.h>
19 
20 namespace gismo {
21 
22 template<typename T>
23 T gsL2Projection<T>::projectGeometry( const gsBasis<T> & basis,
24  const gsGeometry<T> & geometry,
25  gsMatrix<T> & result)
26 {
27  result.clear();
28 
29  gsMultiBasis<T> mb(basis);
30  gsMultiPatch<T> mp;
31  mp.addPatch(geometry);
32 
33  gsExprAssembler<T> A(1,1);
34  A.setIntegrationElements(mb);
35  space u = A.getSpace(mb,mp.targetDim());
36  auto f = A.getCoeff(mp);
37  geometryMap G = A.getMap(mp);
38 
39  u.setup(-1);
40  A.initSystem();
41 
42  // assemble system
43  A.assemble(u*u.tr() * meas(G),u * f * meas(G));
44 
45  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<T>::get( "SimplicialLDLT" );
46  solver->compute(A.matrix());
47  result = solver->solve(A.rhs());
48 
49  solution sol = A.getSolution(u, result);
50  gsExprEvaluator<> ev(A);
51  return ev.integral((sol-f).sqNorm() * meas(G));
52 }
53 
54 
55 template<typename T>
56 T gsL2Projection<T>::projectGeometry( const gsMultiBasis<T> & basis,
57  const gsFunctionSet<T> & geometry,
58  gsMultiPatch<T> & result)
59 {
60  result.clear();
61 
62  gsExprAssembler<T> A(1,1);
63  gsMatrix<T> solVector;
64 
65  A.setIntegrationElements(basis);
66  space u = A.getSpace(basis,geometry.targetDim());
67  solution sol = A.getSolution(u, solVector);
68  auto f = A.getCoeff(geometry);
69  geometryMap G = A.getMap(geometry);
70 
71  u.setup(0);
72  A.initSystem();
73 
74  // assemble system
75  A.assemble(u*u.tr() * meas(G),u * f * meas(G));
76 
77  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<real_t>::get( "SimplicialLDLT" );
78  solver->compute(A.matrix());
79  solVector = solver->solve(A.rhs());
80 
81  gsExprEvaluator<> ev(A);
82  sol.extract(result);
83  result.computeTopology();
84  result.closeGaps();
85 
86  return ev.integral((sol-f).sqNorm() * meas(G));
87 }
88 
89 template<typename T>
90 T gsL2Projection<T>::projectGeometry( const gsMultiBasis<T> & basis,
91  const gsFunctionSet<T> & geometry,
92  gsMatrix<T> & result)
93 {
94  gsExprAssembler<T> A(1,1);
95  A.setIntegrationElements(basis);
96  space u = A.getSpace(basis,geometry.targetDim());
97  auto f = A.getCoeff(geometry);
98  geometryMap G = A.getMap(geometry);
99 
100  u.setup(-1);
101  A.initSystem();
102 
103  // assemble system
104  A.assemble(u*u.tr() * meas(G),u * f * meas(G));
105 
106  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<T>::get( "SimplicialLDLT" );
107  solver->compute(A.matrix());
108  result = solver->solve(A.rhs());
109 
110  solution sol = A.getSolution(u, result);
111  gsExprEvaluator<> ev(A);
112  return ev.integral((sol-f).sqNorm() * meas(G));
113 }
114 
115 template<typename T>
116 T gsL2Projection<T>::projectGeometry( const gsMultiBasis<T> & intbasis,
117  const gsMappedBasis<2,T> & basis,
118  const gsFunctionSet<T> & geometry,
119  gsMatrix<T> & result)
120 {
121  gsExprAssembler<T> A(1,1);
122 
123  A.setIntegrationElements(intbasis);
124  space u = A.getSpace(basis,geometry.targetDim());
125  auto f = A.getCoeff(geometry);
126  geometryMap G = A.getMap(geometry);
127 
128  u.setup(-1);
129  A.initSystem();
130 
131  // assemble system
132  A.assemble(u*u.tr()*meas(G),u * f*meas(G));
133 
134  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<real_t>::get( "SimplicialLDLT" );
135  solver->compute(A.matrix());
136  result = solver->solve(A.rhs());
137 
138  solution sol = A.getSolution(u, result);
139  gsExprEvaluator<> ev(A);
140  return ev.integral((sol-f).sqNorm() * meas(G));
141 }
142 
143 template<typename T>
144 T gsL2Projection<T>::projectFunction( const gsMultiBasis<T> & basis,
145  const gsFunctionSet<T> & source,
146  const gsMultiPatch<T> & geometry,
147  gsMultiPatch<T> & result)
148 {
149  result.clear();
150 
151  gsExprAssembler<T> A(1,1);
152  gsMatrix<T> solVector;
153 
154  A.setIntegrationElements(basis);
155  space u = A.getSpace(basis,source.targetDim());
156  auto f = A.getCoeff(source);
157  solution sol = A.getSolution(u, solVector);
158  geometryMap G = A.getMap(geometry);
159 
160  u.setup(-1);
161  A.initSystem();
162 
163  // assemble system
164  A.assemble(u*u.tr() * meas(G),u * f * meas(G));
165 
166  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<real_t>::get( "SimplicialLDLT" );
167  solver->compute(A.matrix());
168  solVector = solver->solve(A.rhs());
169 
170  sol.extract(result);
171  result.computeTopology();
172  result.closeGaps();
173  gsExprEvaluator<> ev(A);
174  return ev.integral((sol-f).sqNorm() * meas(G));
175 }
176 
177 template<typename T>
178 T gsL2Projection<T>::projectFunction( const gsMultiBasis<T> & basis,
179  const gsFunctionSet<T> & source,
180  const gsMultiPatch<T> & geometry,
181  gsMatrix<T> & result)
182 {
183  gsExprAssembler<T> A(1,1);
184 
185  A.setIntegrationElements(basis);
186  space u = A.getSpace(basis,source.targetDim());
187  auto f = A.getCoeff(source);
188  solution sol = A.getSolution(u, result);
189  geometryMap G = A.getMap(geometry);
190 
191  u.setup(-1);
192  A.initSystem();
193 
194  // assemble system
195  A.assemble(u*u.tr() * meas(G),u * f * meas(G));
196 
197  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<real_t>::get( "SimplicialLDLT" );
198  solver->compute(A.matrix());
199  result = solver->solve(A.rhs());
200  gsExprEvaluator<> ev(A);
201  return ev.integral((sol-f).sqNorm() * meas(G));
202 }
203 
204 
205 template<typename T>
206 T gsL2Projection<T>::projectFunction( const gsMultiBasis<T> & intbasis,
207  const gsMappedBasis<2,T>& basis,
208  const gsFunctionSet<T> & source,
209  const gsMultiPatch<T> & geometry,
210  gsMatrix<T> & result)
211 {
212  gsExprAssembler<T> A(1,1);
213 
214  A.setIntegrationElements(intbasis);
215  space u = A.getSpace(basis,source.targetDim());
216  auto f = A.getCoeff(source);
217  geometryMap G = A.getMap(geometry);
218 
219  u.setup(-1);
220  A.initSystem();
221 
222  // assemble system
223  A.assemble(u*u.tr()*meas(G),u * f *meas(G));
224 
225  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<real_t>::get( "SimplicialLDLT" );
226  solver->compute(A.matrix());
227  result = solver->solve(A.rhs());
228 
229  solution sol = A.getSolution(u, result);
230  gsExprEvaluator<> ev(A);
231  return ev.integral((sol-f).sqNorm() * meas(G));
232 }
233 
234 template<typename T>
235 T gsL2Projection<T>::projectGeometryBoundaries(const gsMultiBasis<T> & basis,
236  const gsMultiPatch<T> & geometry,
237  gsMultiPatch<T> & result)
238 {
239  result.clear();
240 
241  gsExprAssembler<T> A(1,1);
242  gsMatrix<T> solVector;
243 
244  A.setIntegrationElements(basis);
245  space u = A.getSpace(basis,geometry.geoDim());
246  solution sol = A.getSolution(u, solVector);
247  geometryMap G = A.getMap(geometry);
248 
249  std::vector<gsMultiPatch<T>> coords(geometry.geoDim());
250  for (index_t p=0; p!=geometry.geoDim(); p++)
251  coords[p] = geometry.coord(p);
252 
253  gsBoundaryConditions<T> bc;
254  bc.setGeoMap(geometry);
255  for (index_t p=0; p!=geometry.geoDim(); p++)
256  for (typename gsMultiPatch<T>::const_biterator bit = geometry.bBegin(); bit != geometry.bEnd(); ++bit)
257  bc.addCondition(*bit, condition_type::dirichlet,static_cast<gsFunctionSet<T>*>(&coords[p]), 0, true, p);
258 
259  u.setup(bc, dirichlet::l2Projection, 0);
260 
261  A.initSystem();
262 
263  // assemble system
264  A.assemble(u*u.tr()*meas(G),u * G*meas(G));
265 
266  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<real_t>::get( "SimplicialLDLT" );
267  solver->compute(A.matrix());
268  solVector = solver->solve(A.rhs());
269 
270  sol.extract(result);
271  result.computeTopology();
272  result.closeGaps();
273  gsExprEvaluator<> ev(A);
274  return ev.integral((sol-G).sqNorm() * meas(G));
275 }
276 
277 
278 template<typename T>
279 T gsL2Projection<T>::projectGeometryPenalty(const gsMultiBasis<T> & basis,
280  const gsMultiPatch<T> & geometry,
281  gsMultiPatch<T> & result,
282  T penalty)
283 {
284  result.clear();
285 
286  gsExprAssembler<T> A(1,1);
287  gsMatrix<T> solVector;
288 
289  A.setIntegrationElements(basis);
290  space u = A.getSpace(basis,geometry.geoDim());
291  solution sol = A.getSolution(u, solVector);
292  geometryMap G = A.getMap(geometry);
293  auto Gvar = A.getCoeff(geometry,G);
294  element el = A.getElement();
295 
296  std::vector<gsMultiPatch<T>> coords(geometry.geoDim());
297  for (index_t p=0; p!=geometry.geoDim(); p++)
298  coords[p] = geometry.coord(p);
299 
300  gsBoundaryConditions<T> bc;
301  bc.setGeoMap(geometry);
302  for (index_t p=0; p!=geometry.geoDim(); p++)
303  for (typename gsMultiPatch<T>::const_biterator bit = geometry.bBegin(); bit != geometry.bEnd(); ++bit)
304  bc.addCondition(*bit, condition_type::dirichlet,static_cast<gsFunctionSet<T>*>(&coords[p]), 0, true, p);
305 
306 
307  u.setup(bc, dirichlet::l2Projection, -1);
308 
309 
310  std::vector<boundaryInterface> iFace;
311  for (typename gsMultiPatch<T>::const_iiterator iit = geometry.iBegin(); iit != geometry.iEnd(); ++iit)
312  iFace.push_back(*iit);
313 
314  A.initSystem();
315 
316  // assemble system
317  A.assemble(u*u.tr()*meas(G),u * Gvar*meas(G));
318  A.assembleBdr(
319  bc.get("Weak Dirichlet"),
320  -(penalty / el.area(G) * u * u.tr()) * tv(G).norm()
321  );
322 
323  A.assembleIfc(iFace,
324  penalty / el.area(G) * u.left() * u.left().tr() * tv(G).norm()
325  ,
326  - penalty / el.area(G) * u.right()* u.left() .tr() * tv(G).norm()
327  ,
328  - penalty / el.area(G) * u.left() * u.right().tr() * tv(G).norm()
329  ,
330  penalty / el.area(G) * u.right()* u.right().tr() * tv(G).norm()
331  );
332 
333  typename gsSparseSolver<T>::uPtr solver = gsSparseSolver<real_t>::get( "SimplicialLDLT" );
334  solver->compute(A.matrix());
335  solVector = solver->solve(A.rhs());
336 
337  sol.extract(result);
338  result.computeTopology();
339  result.closeGaps();
340  gsExprEvaluator<> ev(A);
341  return ev.integral((Gvar-G).sqNorm() * meas(G));
342 }
343 
344 
345 
346 } // gismo
gismo::gsMultiPatch::addPatch
index_t addPatch(typename gsGeometry< T >::uPtr g)
Add a patch from a gsGeometry&lt;T&gt;::uPtr.
Definition: gsMultiPatch.hpp:210
gismo::gsGeometry
Abstract base class representing a geometry map.
Definition: gsGeometry.h:92
gismo::gsExprEvaluator::integral
T integral(const expr::_expr< E > &expr)
Calculates the integral of the expression expr on the whole integration domain.
Definition: gsExprEvaluator.h:154
gismo::gsExprAssembler
Definition: gsExprAssembler.h:30
gismo::gsMultiPatch::clear
void clear()
Clear (delete) all patches.
Definition: gsMultiPatch.h:319
gismo::condition_type::dirichlet
Dirichlet type.
Definition: gsBoundaryConditions.h:31
gismo::gsMultiPatch::geoDim
short_t geoDim() const
Dimension of the geometry (must match for all patches).
Definition: gsMultiPatch.hpp:149
gsBoundary.h
Provides structs and classes related to interfaces and boundaries.
gismo::gsBoxTopology::iBegin
const_iiterator iBegin() const
Definition: gsBoxTopology.h:119
gismo::gsL2Projection::projectGeometryPenalty
static T projectGeometryPenalty(const gsMultiBasis< T > &basis, const gsMultiPatch< T > &geometry, gsMultiPatch< T > &result, T penalty=1e3)
Projects a source geometry onto basis and returns it in result. Penalizes interfaces and boundaries...
Definition: gsL2Projection.hpp:279
gismo::gsBoxTopology::bBegin
const_biterator bBegin() const
Definition: gsBoxTopology.h:139
gsSparseSolver.h
abstract interfaces for solvers and wrapper around Eigen solvers
gismo::expr::tv
EIGEN_STRONG_INLINE tangent_expr< T > tv(const gsGeometryMap< T > &u)
The tangent boundary vector of a geometry map in 2D.
Definition: gsExpressions.h:4515
gismo::gsExprAssembler::assemble
void assemble(const expr &...args)
Adds the expressions args to the system matrix/rhs The arguments are considered as integrals over the...
Definition: gsExprAssembler.h:756
index_t
#define index_t
Definition: gsConfig.h:32
gismo::gsL2Projection::projectGeometry
static T projectGeometry(const gsBasis< T > &basis, const gsGeometry< T > &geometry, gsMatrix< T > &result)
Projects a source geometry onto basis and returns it in result.
Definition: gsL2Projection.hpp:23
gismo::gsBoundaryConditions::get
bcRefList get(const std::string &label, const short_t unk=0, int comp=-1) const
Definition: gsBoundaryConditions.h:420
gismo::gsFunctionSet::targetDim
virtual short_t targetDim() const
Dimension of the target space.
Definition: gsFunctionSet.h:560
gismo::gsExprAssembler::rhs
const gsMatrix< T > & rhs() const
Returns the right-hand side vector(s)
Definition: gsExprAssembler.h:129
gsExprEvaluator.h
Generic expressions evaluator.
gismo::gsExprAssembler::matrix
const gsSparseMatrix< T > & matrix() const
Returns the left-hand global matrix.
Definition: gsExprAssembler.h:116
gismo::gsL2Projection::projectGeometryBoundaries
static T projectGeometryBoundaries(const gsMultiBasis< T > &basis, const gsMultiPatch< T > &geometry, gsMultiPatch< T > &result)
Projects a source geometry onto basis and returns it in result. Fixes the boundaries.
Definition: gsL2Projection.hpp:235
gismo::gsMultiBasis
Holds a set of patch-wise bases and their topology information.
Definition: gsMultiBasis.h:36
gismo::gsMultiPatch::targetDim
short_t targetDim() const
Dimension of the target space.
Definition: gsMultiPatch.h:192
gismo::expr::gsFeSpace
Definition: gsDirichletValues.h:23
gismo::gsFunctionSet
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition: gsFuncData.h:23
gismo::gsBoxTopology::iEnd
const_iiterator iEnd() const
Definition: gsBoxTopology.h:124
gismo::gsExprAssembler::getMap
geometryMap getMap(const gsFunctionSet< T > &g)
Registers g as an isogeometric geometry map and return a handle to it.
Definition: gsExprAssembler.h:161
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition: gsMultiPatch.h:33
gismo::gsExprEvaluator
Generic evaluator of isogeometric expressions.
Definition: gsExprEvaluator.h:38
gismo::gsExprAssembler::getSpace
space getSpace(const gsFunctionSet< T > &mp, index_t dim=1, index_t id=0)
Definition: gsExprAssembler.h:166
gismo::gsMultiPatch::closeGaps
void closeGaps(T tol=1e-4)
Attempt to close gaps between the interfaces. Assumes that the topology is computed, ie. computeTopology() has been called.
Definition: gsMultiPatch.hpp:566
gismo::gsBoundaryConditions
Class containing a set of boundary conditions.
Definition: gsBoundaryConditions.h:341
gismo::gsBoundaryConditions::addCondition
void addCondition(int p, boxSide s, condition_type::type t, gsFunctionSet< T > *f, short_t unknown=0, bool parametric=false, int comp=-1)
Adds another boundary condition.
Definition: gsBoundaryConditions.h:650
gsBoundaryConditions.h
Provides gsBoundaryConditions class.
gismo::gsSparseSolver
Abstract class for solvers. The solver interface is base on 3 methods: -compute set the system matrix...
Definition: gsSparseSolver.h:66
gismo::expr::meas
EIGEN_STRONG_INLINE meas_expr< T > meas(const gsGeometryMap< T > &G)
The measure of a geometry map.
Definition: gsExpressions.h:4557
gismo::gsBoundaryConditions::setGeoMap
void setGeoMap(const gsFunctionSet< T > &gm)
Set the geometry map to evaluate boundary conditions.
Definition: gsBoundaryConditions.h:916
gismo::gsExprAssembler::setIntegrationElements
void setIntegrationElements(const gsMultiBasis< T > &mesh)
Sets the domain of integration.
Definition: gsExprAssembler.h:136
gismo::gsBoxTopology::bEnd
const_biterator bEnd() const
Definition: gsBoxTopology.h:144
gismo::gsExprAssembler::getSolution
solution getSolution(const expr::gsFeSpace< T > &s, gsMatrix< T > &cf) const
Registers a representation of a solution variable from space s, based on the vector cf...
Definition: gsExprAssembler.h:274
gismo::gsExprAssembler::initSystem
void initSystem(const index_t numRhs=1)
Initializes the sparse system (sparse matrix and rhs)
Definition: gsExprAssembler.h:290
gismo::gsExprAssembler::getCoeff
variable getCoeff(const gsFunctionSet< T > &func)
Definition: gsExprAssembler.h:260
gsExprAssembler.h
Generic expressions matrix assembly.
gismo::gsExprAssembler::assembleBdr
void assembleBdr(const bcRefList &BCs, expr &...args)
Adds the expressions args to the system matrix/rhs The arguments are considered as integrals over the...
Definition: gsExprAssembler.h:837
gismo::gsL2Projection::projectFunction
static T projectFunction(const gsMultiBasis< T > &basis, const gsFunctionSet< T > &source, const gsMultiPatch< T > &geometry, gsMatrix< T > &result)
Projects a function on a basis.
Definition: gsL2Projection.hpp:178
gismo::gsBasis
A basis represents a family of scalar basis functions defined over a common parameter domain...
Definition: gsBasis.h:78
gismo::gsMultiPatch::computeTopology
bool computeTopology(T tol=1e-4, bool cornersOnly=false, bool tjunctions=false)
Attempt to compute interfaces and boundaries automatically.
Definition: gsMultiPatch.hpp:366