Line data Source code
1 : /*
2 : * Copyright (c) 2020: G-CSC, Goethe University Frankfurt
3 : * Author: Dmitry Logashenko
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : /*
34 : * Global assembling of the problems with the embedded boundary.
35 : */
36 : #ifndef __H__UG__PLUGINS__D3F__EMBASS__
37 : #define __H__UG__PLUGINS__D3F__EMBASS__
38 :
39 : #include <vector>
40 :
41 : // ug4 headers
42 : #include "common/common.h"
43 : #include "common/util/smart_pointer.h"
44 : #include "lib_disc/domain_traits.h"
45 : #include "lib_disc/spatial_disc/domain_disc.h"
46 : #include "lib_disc/operator/linear_operator/std_injection.h"
47 : #include "bridge/util_algebra_dependent.h"
48 : #ifdef UG_FOR_LUA
49 : #include "bindings/lua/lua_user_data.h"
50 : #endif
51 :
52 : namespace ug {
53 :
54 : /// Base class for the extrapolation over an embedded boundary
55 : /**
56 : * This class provides an interface for the access to the extrapolation over
57 : * an embedded boundary.
58 : *
59 : * \tparam TDomain domain type
60 : * \tparam TAlgebra algebra type for the functions to extrapolate
61 : */
62 : template <typename TDomain, typename TAlgebra>
63 : class IInterfaceExtrapolation
64 : {
65 : public:
66 :
67 : /// domain type
68 : typedef TDomain domain_type;
69 :
70 : /// algebra type for the functions to extrapolate
71 : typedef TAlgebra algebra_type;
72 :
73 : /// vector type (for the functions to extrapolate)
74 : typedef typename algebra_type::vector_type vector_type;
75 :
76 : /// matrix type
77 : typedef typename algebra_type::matrix_type matrix_type;
78 :
79 : /// dimensionality (the World dimension)
80 : static const int dim = domain_type::dim;
81 :
82 : /// Constructor
83 : IInterfaceExtrapolation () {}
84 :
85 : /// Destructor
86 0 : virtual ~IInterfaceExtrapolation () {}
87 :
88 : /// checks whether the element is intersected by the interface, or what, and prepares the data
89 : virtual int check_elem_lsf
90 : (
91 : size_t n_co, ///< number of the corners of the element
92 : GridObject * pElem, ///< the element to process
93 : int si, ///< subset of the element
94 : int g_level, ///< grid level of the element
95 : bool use_hanging, ///< if there can be hanging nodes
96 : const MathVector<dim> vCornerCoords [], ///< coordinates of the corners of the element
97 : number time ///< the phisical time
98 : )
99 : {UG_THROW ("IInterfaceExtrapolation: Virtual functions are not implemented in the base class.");}
100 :
101 : /// (slower version) checks whether the element is intersected by the interface, or what, and prepares the data
102 : virtual int check_elem_lsf
103 : (
104 : size_t n_co, ///< number of the corners of the element
105 : GridObject * pElem, ///< the element to process
106 : int si, ///< subset of the element
107 : bool use_hanging, ///< if there can be hanging nodes
108 : const MathVector<dim> vCornerCoords [], ///< coordinates of the corners of the element
109 : number time ///< the phisical time
110 : )
111 : {UG_THROW ("IInterfaceExtrapolation: Virtual functions are not implemented in the base class.");}
112 :
113 : /// extrapolates a component of the solution to the vertices behind the interface (w.r.t. a base corner)
114 : virtual void extrapolate_by_lsf
115 : (
116 : size_t num_co, ///< number of the corners
117 : size_t base_co, ///< the base corner
118 : number * u, ///< nodal values to extrapolate
119 : size_t fct ///< index of the function (to identify to type of the extrapolation)
120 : ) const
121 : {UG_THROW ("IInterfaceExtrapolation: Virtual functions are not implemented in the base class.");}
122 :
123 : /// extrapolates a component of the solution to the vertices behind the interface (by averaging)
124 : virtual void extrapolate_by_lsf
125 : (
126 : size_t num_co, ///< number of the corners
127 : number * u, ///< nodal values to extrapolate
128 : size_t fct ///< index of the function (to identify to type of the extrapolation)
129 : ) const
130 : {UG_THROW ("IInterfaceExtrapolation: Virtual functions are not implemented in the base class.");}
131 :
132 : /// returns true if the corner is "inside" (use after check_elem_lsf)
133 : virtual bool corner_inside
134 : (
135 : size_t co ///< the corner
136 : ) const
137 : {UG_THROW ("IInterfaceExtrapolation: Virtual functions are not implemented in the base class.");}
138 :
139 : /// returns the effective value of the LSF at a corner (use after check_elem_lsf)
140 : virtual number lsf_at
141 : (
142 : size_t co ///< the corner
143 : ) const
144 : {UG_THROW ("IInterfaceExtrapolation: Virtual functions are not implemented in the base class.");}
145 :
146 : }; // class IInterfaceExtrapolation
147 :
148 : /// Global assembler based on the ghost-fluid method with a level-set function
149 : /**
150 : * Template class of the global assembler based on the ghost-fluid method
151 : * with a piecewise linear level-set function.
152 : *
153 : * \tparam TDomain domain type
154 : * \tparam TAlgebra algebra type
155 : * \tparam TExtrapolation extrapolation class
156 : */
157 : template <typename TDomain, typename TAlgebra, typename TExtrapolation>
158 : class LSGFGlobAssembler
159 : {
160 : public:
161 :
162 : /// Domain type
163 : typedef TDomain domain_type;
164 :
165 : /// Algebra type
166 : typedef TAlgebra algebra_type;
167 :
168 : /// type of approximation space
169 : typedef ApproximationSpace<domain_type> approx_space_type;
170 :
171 : /// Vector type in the algebra
172 : typedef typename algebra_type::vector_type vector_type;
173 :
174 : /// Matrix type in the algebra
175 : typedef typename algebra_type::matrix_type matrix_type;
176 :
177 : /// Extrapolation type
178 : typedef TExtrapolation extrapolation_type;
179 :
180 : /// Grid function type for the LSF
181 : typedef typename extrapolation_type::ls_grid_func_type ls_grid_func_type;
182 :
183 : /// world dimension
184 : static const int dim = TDomain::dim;
185 :
186 : ////////////////////////////////////////////////////////////////////////////////
187 : // Constructor/destructor
188 : ////////////////////////////////////////////////////////////////////////////////
189 :
190 : public:
191 :
192 : /// class constructor (may not have any arguments!)
193 : LSGFGlobAssembler () : m_bAssembleOnlyCut(false) {};
194 :
195 : /// virtual destructor
196 : virtual ~LSGFGlobAssembler () {};
197 :
198 : ////////////////////////////////////////////////////////////////////////////////
199 : // Assembling tools
200 : ////////////////////////////////////////////////////////////////////////////////
201 :
202 : public:
203 :
204 : template <typename TElem, typename TIterator>
205 : void
206 : AssembleStiffnessMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
207 : ConstSmartPtr<domain_type> spDomain,
208 : ConstSmartPtr<DoFDistribution> dd,
209 : TIterator iterBegin,
210 : TIterator iterEnd,
211 : int si, bool bNonRegularGrid,
212 : matrix_type& A,
213 : const vector_type& u,
214 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
215 :
216 : template <typename TElem, typename TIterator>
217 : void
218 : AssembleMassMatrix( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
219 : ConstSmartPtr<domain_type> spDomain,
220 : ConstSmartPtr<DoFDistribution> dd,
221 : TIterator iterBegin,
222 : TIterator iterEnd,
223 : int si, bool bNonRegularGrid,
224 : matrix_type& M,
225 : const vector_type& u,
226 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
227 :
228 : template <typename TElem, typename TIterator>
229 : void
230 : AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
231 : ConstSmartPtr<domain_type> spDomain,
232 : ConstSmartPtr<DoFDistribution> dd,
233 : TIterator iterBegin,
234 : TIterator iterEnd,
235 : int si, bool bNonRegularGrid,
236 : matrix_type& J,
237 : const vector_type& u,
238 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
239 :
240 : template <typename TElem, typename TIterator>
241 : void
242 : AssembleJacobian( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
243 : ConstSmartPtr<domain_type> spDomain,
244 : ConstSmartPtr<DoFDistribution> dd,
245 : TIterator iterBegin,
246 : TIterator iterEnd,
247 : int si, bool bNonRegularGrid,
248 : matrix_type& J,
249 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
250 : number s_a0,
251 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
252 :
253 : template <typename TElem, typename TIterator>
254 : void
255 : AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
256 : ConstSmartPtr<domain_type> spDomain,
257 : ConstSmartPtr<DoFDistribution> dd,
258 : TIterator iterBegin,
259 : TIterator iterEnd,
260 : int si, bool bNonRegularGrid,
261 : vector_type& d,
262 : const vector_type& u,
263 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
264 :
265 : template <typename TElem, typename TIterator>
266 : void
267 : AssembleDefect( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
268 : ConstSmartPtr<domain_type> spDomain,
269 : ConstSmartPtr<DoFDistribution> dd,
270 : TIterator iterBegin,
271 : TIterator iterEnd,
272 : int si, bool bNonRegularGrid,
273 : vector_type& d,
274 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
275 : const std::vector<number>& vScaleMass,
276 : const std::vector<number>& vScaleStiff,
277 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
278 :
279 : template <typename TElem, typename TIterator>
280 : void
281 : AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
282 : ConstSmartPtr<domain_type> spDomain,
283 : ConstSmartPtr<DoFDistribution> dd,
284 : TIterator iterBegin,
285 : TIterator iterEnd,
286 : int si, bool bNonRegularGrid,
287 : matrix_type& A,
288 : vector_type& rhs,
289 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
290 :
291 : template <typename TElem, typename TIterator>
292 : void
293 : AssembleLinear( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
294 : ConstSmartPtr<domain_type> spDomain,
295 : ConstSmartPtr<DoFDistribution> dd,
296 : TIterator iterBegin,
297 : TIterator iterEnd,
298 : int si, bool bNonRegularGrid,
299 : matrix_type& A,
300 : vector_type& rhs,
301 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
302 : const std::vector<number>& vScaleMass,
303 : const std::vector<number>& vScaleStiff,
304 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
305 :
306 : ////////////////////////////////////////////////////////////////////////////////
307 : // Assemble Rhs: it cannot be done for the ghost-fluid method independently of the matrix
308 : ////////////////////////////////////////////////////////////////////////////////
309 :
310 : public:
311 :
312 : template <typename TElem, typename TIterator>
313 : static void
314 : AssembleRhs( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
315 : ConstSmartPtr<domain_type> spDomain,
316 : ConstSmartPtr<DoFDistribution> dd,
317 : TIterator iterBegin,
318 : TIterator iterEnd,
319 : int si, bool bNonRegularGrid,
320 : vector_type& rhs,
321 : const vector_type& u,
322 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
323 : {
324 : UG_THROW ("LSGFGlobAssembler::AssembleRhs: Cannot assemble the RHS in GF independently of the matrix");
325 : }
326 :
327 : template <typename TElem, typename TIterator>
328 : static void
329 : AssembleRhs( const std::vector<IElemDisc<domain_type>*>& vElemDisc,
330 : ConstSmartPtr<domain_type> spDomain,
331 : ConstSmartPtr<DoFDistribution> dd,
332 : TIterator iterBegin,
333 : TIterator iterEnd,
334 : int si, bool bNonRegularGrid,
335 : vector_type& rhs,
336 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
337 : const std::vector<number>& vScaleMass,
338 : const std::vector<number>& vScaleStiff,
339 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner)
340 : {
341 : UG_THROW ("LSGFGlobAssembler::AssembleRhs: Cannot assemble the RHS in GF independently of the matrix");
342 : }
343 :
344 : ////////////////////////////////////////////////////////////////////////////////
345 : // Prepare and Finish Timestep: these version merely skip the outer elements
346 : ////////////////////////////////////////////////////////////////////////////////
347 :
348 : public:
349 : /**
350 : * This function prepares the global discretization for a time-stepping scheme
351 : * by calling the "prepare_timestep" methods of all passed element
352 : * discretizations.
353 : *
354 : * \param[in] vElemDisc element discretizations
355 : * \param[in] dd DoF Distribution
356 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
357 : * \param[in] vSol current and previous solutions
358 : * \param[in] spAssTuner assemble adapter
359 : */
360 : void PrepareTimestep
361 : (
362 : const std::vector<IElemDisc<domain_type>*>& vElemDisc,
363 : ConstSmartPtr<DoFDistribution> dd,
364 : bool bNonRegularGrid,
365 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
366 : number future_time,
367 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner
368 : );
369 :
370 : template <typename TElem, typename TIterator>
371 : void
372 : PrepareTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
373 : ConstSmartPtr<domain_type> spDomain,
374 : ConstSmartPtr<DoFDistribution> dd,
375 : TIterator iterBegin,
376 : TIterator iterEnd,
377 : int si, bool bNonRegularGrid,
378 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
379 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
380 : /**
381 : * This function finishes the global discretization for a time-stepping scheme
382 : * by calling the "finish_timestep" methods of all passed element
383 : * discretizations.
384 : *
385 : * \param[in] vElemDisc element discretizations
386 : * \param[in] dd DoF Distribution
387 : * \param[in] bNonRegularGrid flag to indicate if non regular grid is used
388 : * \param[in] vSol current and previous solutions
389 : * \param[in] spAssTuner assemble adapter
390 : */
391 : void FinishTimestep
392 : (
393 : const std::vector<IElemDisc<domain_type>*>& vElemDisc,
394 : ConstSmartPtr<DoFDistribution> dd,
395 : bool bNonRegularGrid,
396 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
397 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner
398 : );
399 :
400 : template <typename TElem, typename TIterator>
401 : void
402 : FinishTimestepElem(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
403 : ConstSmartPtr<domain_type> spDomain,
404 : ConstSmartPtr<DoFDistribution> dd,
405 : TIterator iterBegin,
406 : TIterator iterEnd,
407 : int si, bool bNonRegularGrid,
408 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
409 : ConstSmartPtr<AssemblingTuner<TAlgebra> > spAssTuner);
410 :
411 : template <typename TElem, typename TIterator>
412 : void
413 : InitAllExports(const std::vector<IElemDisc<domain_type>*>& vElemDisc,
414 : ConstSmartPtr<DoFDistribution> dd,
415 : TIterator iterBegin,
416 : TIterator iterEnd,
417 : int si, bool bNonRegularGrid, bool bAsTimeDependent);
418 :
419 : ////////////////////////////////////////////////////////////////////////////////
420 : // Error estimators: Not implemented for the ghost-fluid method
421 : ////////////////////////////////////////////////////////////////////////////////
422 :
423 : public:
424 :
425 : template <typename TElem, typename TIterator>
426 : static void
427 : AssembleErrorEstimator
428 : (
429 : const std::vector<IElemError<domain_type>*>& vElemDisc,
430 : ConstSmartPtr<domain_type> spDomain,
431 : ConstSmartPtr<DoFDistribution> dd,
432 : TIterator iterBegin,
433 : TIterator iterEnd,
434 : int si,
435 : bool bNonRegularGrid,
436 : const vector_type& u
437 : )
438 : {
439 : UG_THROW ("AssembleErrorEstimator: No error estimator implemented for the Ghost-Fluid method.");
440 : }
441 :
442 : template <typename TElem, typename TIterator>
443 : static void
444 : AssembleErrorEstimator
445 : (
446 : const std::vector<IElemError<domain_type>*>& vElemDisc,
447 : ConstSmartPtr<domain_type> spDomain,
448 : ConstSmartPtr<DoFDistribution> dd,
449 : TIterator iterBegin,
450 : TIterator iterEnd,
451 : int si,
452 : bool bNonRegularGrid,
453 : std::vector<number> vScaleMass,
454 : std::vector<number> vScaleStiff,
455 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol
456 : )
457 : {
458 : UG_THROW ("AssembleErrorEstimator: No error estimator implemented for the Ghost-Fluid method.");
459 : }
460 :
461 : ////////////////////////////////////////////////////////////////////////////////
462 : // Data
463 : ////////////////////////////////////////////////////////////////////////////////
464 :
465 : private:
466 :
467 : extrapolation_type m_extrapol; ///< the extrapolation at the interface
468 :
469 : /**
470 : * The following flag is used only for research purposes (measuring the volume
471 : * sources/sinks at the embedded interface etc.). It switches off assembing
472 : * in all inner elements (i.e. makes the procedures to assemble only in the cut
473 : * elements:
474 : */
475 : bool m_bAssembleOnlyCut;
476 :
477 : public:
478 :
479 : /// set the level-set function and check it
480 : void set_LSF
481 : (
482 : SmartPtr<ls_grid_func_type> spLSF ///< the function to set
483 : )
484 : {
485 : m_extrapol.set_LSF (spLSF);
486 : }
487 :
488 : /// prepares the boundary conditions at the interface: sets all them to Dirichlet-0
489 : void prepare_interface_bc
490 : (
491 : SmartPtr<approx_space_type> spApproxSpace ///< the approximation space of the domain discretization
492 : )
493 : {
494 : m_extrapol.prepare_interface_bc (spApproxSpace);
495 : }
496 :
497 : /// adds a Dirichlet BC with a given value on the interface
498 : void set_Dirichlet_on_if_for
499 : (
500 : SmartPtr<approx_space_type> spApproxSpace, ///< the approximation space of the domain discretization
501 : const char* fct_name, ///< function to impose the condition for
502 : number value ///< the Dirichlet value
503 : )
504 : {
505 : m_extrapol.set_Dirichlet_for (spApproxSpace, fct_name, value);
506 : }
507 :
508 : /// adds a Dirichlet BC with a given value on the interface
509 : void set_Dirichlet_on_if_for
510 : (
511 : SmartPtr<approx_space_type> spApproxSpace, ///< the approximation space of the domain discretization
512 : const char* fct_name, ///< function to impose the condition for
513 : SmartPtr<CplUserData<number, dim> > func ///< the Dirichlet function
514 : )
515 : {
516 : m_extrapol.set_Dirichlet_for (spApproxSpace, fct_name, func);
517 : }
518 :
519 : /// adds a "plain" Dirichlet BC with a given value on the interface
520 : void set_plain_Dirichlet_on_if_for
521 : (
522 : SmartPtr<approx_space_type> spApproxSpace, ///< the approximation space of the domain discretization
523 : const char* fct_name, ///< function to impose the condition for
524 : SmartPtr<CplUserData<number, dim> > func ///< the Dirichlet function
525 : )
526 : {
527 : m_extrapol.set_plain_Dirichlet_for (spApproxSpace, fct_name, func);
528 : }
529 :
530 : /// adds a "plain" Dirichlet BC with a given value on the interface
531 : void set_plain_Dirichlet_on_if_for
532 : (
533 : SmartPtr<approx_space_type> spApproxSpace, ///< the approximation space of the domain discretization
534 : const char* fct_name, ///< function to impose the condition for
535 : number value ///< the Dirichlet value
536 : )
537 : {
538 : m_extrapol.set_plain_Dirichlet_for (spApproxSpace, fct_name, value);
539 : }
540 :
541 : /// adds a Neumann-0 with on the interface
542 : void set_Neumann0_on_if_for
543 : (
544 : SmartPtr<approx_space_type> spApproxSpace, ///< the approximation space of the domain discretization
545 : const char* fct_name ///< function to impose the condition for
546 : )
547 : {
548 : m_extrapol.set_Neumann0_for (spApproxSpace, fct_name);
549 : }
550 :
551 : /// excludes a (boundary) subsets from the extrapolation
552 : void exclude_subsets
553 : (
554 : SmartPtr<approx_space_type> spApproxSpace, ///< the approximation space of the domain discretization
555 : const char* subset_names ///< names of the subsets to exclude
556 : )
557 : {
558 : m_extrapol.exclude_subsets (spApproxSpace, subset_names);
559 : }
560 :
561 : /// set the "assemble only in cut elements" flag
562 : void set_assemble_only_cut (bool b)
563 : {
564 : m_bAssembleOnlyCut = b;
565 : }
566 :
567 : /// project the level-set function to the coarse levels
568 : void project_LSF ()
569 : {
570 : m_extrapol.project_LSF ();
571 : }
572 :
573 : /// returns the extrapolation
574 : extrapolation_type & extrapolation () {return m_extrapol;}
575 :
576 : /// checks whether the element is intersected by the interface, or what, and prepares the data
577 : int check_elem_lsf
578 : (
579 : size_t n_co, ///< number of the corners of the element
580 : GridObject * pElem, ///< the element to process
581 : int si, ///< subset of the element
582 : int g_level, ///< grid level of the element
583 : bool use_hanging, ///< if there can be hanging nodes
584 : const MathVector<dim> vCornerCoords [], ///< coordinates of the corners of the element
585 : number time ///< the phisical time
586 : )
587 : {
588 : return m_extrapol.check_elem_lsf (n_co, pElem, si, g_level, use_hanging, vCornerCoords, time);
589 : }
590 :
591 : /// extrapolates a component the solution to the vertices behind the interface (w.r.t. a base corner)
592 : void extrapolate_by_lsf
593 : (
594 : size_t num_co, ///< number of the corners
595 : size_t base_co, ///< the base corner
596 : number * u, ///< nodal values to extrapolate
597 : size_t fct ///< index of the function (to identify to type of the extrapolation)
598 : ) const
599 : {
600 : m_extrapol.extrapolate_by_lsf (num_co, base_co, u, fct);
601 : }
602 :
603 : /// extrapolates a component the solution to the vertices behind the interface (by averaging)
604 : void extrapolate_by_lsf
605 : (
606 : size_t num_co, ///< number of the corners
607 : number * u, ///< nodal values to extrapolate
608 : size_t fct ///< index of the function (to identify to type of the extrapolation)
609 : ) const
610 : {
611 : m_extrapol.extrapolate_by_lsf (num_co, u, fct);
612 : }
613 :
614 : /// returns true if the corner is "inside" (use after check_elem_lsf)
615 : bool corner_inside
616 : (
617 : size_t co ///< the corner
618 : ) const
619 : {return m_extrapol.corner_inside (co);}
620 :
621 : /// returns the effective value of the LSF at a corner (use after check_elem_lsf)
622 : number lsf_at
623 : (
624 : size_t co ///< the corner
625 : ) const
626 : {return m_extrapol.lsf_at (co);}
627 :
628 : /// sets the values at the outer vertices to 0
629 : virtual void clear_outer_values
630 : (
631 : vector_type & d, ///< the vector where to set
632 : const DoFDistribution * dd ///< dof distribution of the grid function to reset
633 : ) const
634 : {m_extrapol.clear_outer_values (d, dd);}
635 :
636 : /// sets the values at the outer vertices to given values
637 : virtual void set_outer_values
638 : (
639 : vector_type & u, ///< the vector where to set
640 : const DoFDistribution * dd, ///< dof distribution of the grid function to reset
641 : number time ///< the physical time
642 : )
643 : {m_extrapol.set_outer_values (u, dd, time);}
644 :
645 : /// sets the matrices at outer vertices to identity
646 : virtual void set_outer_matrices
647 : (
648 : matrix_type & A, ///< the vector where to set
649 : const DoFDistribution * dd ///< dof distribution of the grid function to reset
650 : ) const
651 : {m_extrapol.set_outer_matrices (A, dd);}
652 :
653 : }; // class LSGFGlobAssembler
654 :
655 : /// a special constraint that sets functions and matrices in the outer subdomain to given values
656 : /**
657 : * This class implements a constraint for the LSF based ghost-fluid method
658 : * to set grid functions and matrices in the outer subdomain to given values.
659 : *
660 : * \tparam TDomain domain type
661 : * \tparam TAlgebra algebra type
662 : * \tparam TExtrapolation extrapolation class
663 : */
664 : template <typename TDomain, typename TAlgebra, typename TExtrapolation>
665 : class LSGFConstraint
666 : : public IDomainConstraint<TDomain, TAlgebra>
667 : {
668 : public:
669 :
670 : /// Domain type
671 : typedef TDomain domain_type;
672 :
673 : /// Algebra type
674 : typedef TAlgebra algebra_type;
675 :
676 : /// Vector type in the algebra
677 : typedef typename algebra_type::vector_type vector_type;
678 :
679 : /// Matrix type in the algebra
680 : typedef typename algebra_type::matrix_type matrix_type;
681 :
682 : /// Extrapolation type
683 : typedef TExtrapolation extrapolation_type;
684 :
685 : /// class constructor
686 : LSGFConstraint
687 : (
688 : extrapolation_type & rExtrapolation ///< the GF extrapolation
689 : )
690 : : m_rExtrapolation (rExtrapolation)
691 : {}
692 :
693 : /// virtual destructor
694 : virtual ~LSGFConstraint () {};
695 :
696 : /// sets a unity row for all conductor indices
697 : void adjust_jacobian
698 : (
699 : matrix_type & J,
700 : const vector_type & u,
701 : ConstSmartPtr<DoFDistribution> dd,
702 : int type,
703 : number time = 0.0,
704 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = SPNULL,
705 : const number s_a0 = 1.0
706 : )
707 : {
708 : m_rExtrapolation.set_outer_matrices (J, dd.get ());
709 : }
710 :
711 : /// sets a zero value in the defect for all conductor indices
712 : void adjust_defect
713 : (
714 : vector_type & d,
715 : const vector_type & u,
716 : ConstSmartPtr<DoFDistribution> dd,
717 : int type,
718 : number time = 0.0,
719 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = SPNULL,
720 : const std::vector<number> * vScaleMass = NULL,
721 : const std::vector<number> * vScaleStiff = NULL
722 : )
723 : {
724 : m_rExtrapolation.clear_outer_values (d, dd.get ());
725 : }
726 :
727 : /// sets the value in the solution for all conductor indices
728 : void adjust_solution
729 : (
730 : vector_type & u,
731 : ConstSmartPtr<DoFDistribution> dd,
732 : int type,
733 : number time = 0.0
734 : )
735 : {
736 : m_rExtrapolation.set_outer_values (u, dd.get (), time);
737 : }
738 :
739 : /// sets unity rows in A and dirichlet values in right-hand side b
740 : void adjust_linear
741 : (
742 : matrix_type & A,
743 : vector_type & b,
744 : ConstSmartPtr<DoFDistribution> dd,
745 : int type,
746 : number time = 0.0
747 : )
748 : { // Note that this function is not really used, so it needs not to be optimal.
749 : m_rExtrapolation.set_outer_matrices (A, dd.get ());
750 : m_rExtrapolation.clear_outer_values (b, dd.get ());
751 : }
752 :
753 : /// sets the dirichlet value in the right-hand side
754 : void adjust_rhs
755 : (
756 : vector_type & b,
757 : const vector_type & u,
758 : ConstSmartPtr<DoFDistribution> dd,
759 : int type,
760 : number time = 0.0
761 : )
762 : {
763 : m_rExtrapolation.clear_outer_values (b, dd.get ());
764 : }
765 :
766 : /// returns the type of the constraints
767 : int type () const {return CT_DIRICHLET;}
768 :
769 : private:
770 :
771 : /// Extrapolation in the GF method
772 : extrapolation_type & m_rExtrapolation;
773 : };
774 :
775 : /// domain discretization for the Level-Set Ghost-Fluid method
776 : /**
777 : * This class template is an implementation of the IDomainDiscretization
778 : * interface for the Level-Set Ghost-Fluid method.
779 : *
780 : * \tparam TDomain domain type
781 : * \tparam TAlgebra algebra type
782 : */
783 : template <typename TDomain, typename TAlgebra, typename TExtrapolation>
784 : class LSGFDomainDiscretization
785 : : public DomainDiscretizationBase<TDomain, TAlgebra, LSGFGlobAssembler<TDomain, TAlgebra, TExtrapolation> >,
786 : public IInterfaceExtrapolation<TDomain, TAlgebra>
787 : {
788 : /// Type of the global assembler
789 : typedef LSGFGlobAssembler<TDomain, TAlgebra, TExtrapolation> gass_type;
790 :
791 : /// Type of the base class
792 : typedef DomainDiscretizationBase<TDomain, TAlgebra, gass_type> base_type;
793 :
794 : /// Type of the constraint
795 : typedef LSGFConstraint<TDomain, TAlgebra, TExtrapolation> ls_constraint_type;
796 :
797 : public:
798 : /// Type of Domain
799 : typedef TDomain domain_type;
800 :
801 : /// Type of algebra
802 : typedef TAlgebra algebra_type;
803 :
804 : /// Type of the grid
805 : typedef typename TDomain::grid_type grid_type;
806 :
807 : /// Type of algebra matrix
808 : typedef typename algebra_type::matrix_type matrix_type;
809 :
810 : /// Type of algebra vector
811 : typedef typename algebra_type::vector_type vector_type;
812 :
813 : /// Type of approximation space
814 : typedef ApproximationSpace<TDomain> approx_space_type;
815 :
816 : /// Type of the LSF grid functions
817 : typedef typename gass_type::ls_grid_func_type ls_grid_func_type;
818 :
819 : /// world dimension
820 : static const int dim = TDomain::dim;
821 :
822 : public:
823 : /// default Constructor
824 : LSGFDomainDiscretization (SmartPtr<approx_space_type> pApproxSpace)
825 : : DomainDiscretizationBase<domain_type, algebra_type, gass_type> (pApproxSpace),
826 : m_spLSFGFConstraint (new ls_constraint_type (this->extrapolation ()))
827 : {
828 : // register the constraint
829 : this->add (SmartPtr<IDomainConstraint<domain_type, algebra_type> > (m_spLSFGFConstraint));
830 : // set the default the boundary condtions at the interface
831 : gass_type::prepare_interface_bc (pApproxSpace);
832 : }
833 :
834 : /// virtual destructor
835 : virtual ~LSGFDomainDiscretization () {};
836 :
837 : /// set the level-set function and check it
838 : void set_LSF
839 : (
840 : SmartPtr<ls_grid_func_type> spLSF ///< the function to set
841 : )
842 : {
843 : gass_type::set_LSF (spLSF);
844 : }
845 :
846 : /// sets the Dirichlet boundary condition at the interface for a component of the solution
847 : void set_Dirichlet_on_if_for
848 : (
849 : const char* fct_name, ///< function to impose the condition for
850 : number value ///< the Dirichlet value
851 : )
852 : {
853 : gass_type::set_Dirichlet_on_if_for (base_type::m_spApproxSpace, fct_name, value);
854 : }
855 :
856 : /// sets the Dirichlet boundary condition at the interface for a component of the solution
857 : void set_Dirichlet_on_if_for
858 : (
859 : const char* fct_name, ///< function to impose the condition for
860 : SmartPtr<CplUserData<number, dim> > func ///< the Dirichlet value (as a function)
861 : )
862 : {
863 : gass_type::set_Dirichlet_on_if_for (base_type::m_spApproxSpace, fct_name, func);
864 : }
865 :
866 : /// sets the 'plain' Dirichlet boundary condition at the interface for a component of the solution
867 : void set_plain_Dirichlet_on_if_for
868 : (
869 : const char* fct_name, ///< function to impose the condition for
870 : number value ///< the Dirichlet value
871 : )
872 : {
873 : gass_type::set_plain_Dirichlet_on_if_for (base_type::m_spApproxSpace, fct_name, value);
874 : }
875 :
876 : /// sets the 'plain' Dirichlet boundary condition at the interface for a component of the solution
877 : void set_plain_Dirichlet_on_if_for
878 : (
879 : const char* fct_name, ///< function to impose the condition for
880 : SmartPtr<CplUserData<number, dim> > func ///< the Dirichlet value (as a function)
881 : )
882 : {
883 : gass_type::set_plain_Dirichlet_on_if_for (base_type::m_spApproxSpace, fct_name, func);
884 : }
885 :
886 : #ifdef UG_FOR_LUA
887 :
888 : /// sets the Dirichlet boundary condition at the interface for a component of the solution
889 : void set_Dirichlet_on_if_for
890 : (
891 : const char* fct_name, ///< function to impose the condition for
892 : const char* func_name ///< the Dirichlet value (as a name of the LUA function)
893 : )
894 : {
895 : set_Dirichlet_on_if_for (fct_name, LuaUserDataFactory<number,dim>::create(func_name));
896 : }
897 :
898 : /// sets the Dirichlet boundary condition at the interface for a component of the solution
899 : void set_Dirichlet_on_if_for
900 : (
901 : const char* fct_name, ///< function to impose the condition for
902 : LuaFunctionHandle func ///< the Dirichlet value (as a LUA function)
903 : )
904 : {
905 : set_Dirichlet_on_if_for (fct_name, make_sp(new LuaUserData<number,dim>(func)));
906 : }
907 :
908 : /// sets the Dirichlet boundary condition at the interface for a component of the solution
909 : void set_plain_Dirichlet_on_if_for
910 : (
911 : const char* fct_name, ///< function to impose the condition for
912 : const char* func_name ///< the Dirichlet value (as a name of the LUA function)
913 : )
914 : {
915 : set_plain_Dirichlet_on_if_for (fct_name, LuaUserDataFactory<number,dim>::create(func_name));
916 : }
917 :
918 : /// sets the Dirichlet boundary condition at the interface for a component of the solution
919 : void set_plain_Dirichlet_on_if_for
920 : (
921 : const char* fct_name, ///< function to impose the condition for
922 : LuaFunctionHandle func ///< the Dirichlet value (as a LUA function)
923 : )
924 : {
925 : set_plain_Dirichlet_on_if_for (fct_name, make_sp(new LuaUserData<number,dim>(func)));
926 : }
927 :
928 : #endif
929 :
930 : /// sets the Neumann-0 boundary condition at the interface for a component of the solution
931 : void set_Neumann0_on_if_for
932 : (
933 : const char* fct_name ///< function to impose the condition for
934 : )
935 : {
936 : gass_type::set_Neumann0_on_if_for (base_type::m_spApproxSpace, fct_name);
937 : }
938 :
939 : /// excludes a (boundary) subsets from the extrapolation
940 : void exclude_subsets
941 : (
942 : const char* subset_names ///< names of the subsets to exclude
943 : )
944 : {
945 : gass_type::exclude_subsets (base_type::m_spApproxSpace, subset_names);
946 : }
947 :
948 : /// set the "assemble only in cut elements" flag
949 : void set_assemble_only_cut (bool b)
950 : {
951 : gass_type::set_assemble_only_cut (b);
952 : }
953 :
954 : /// project the level-set function to the coarse levels
955 : void project_LSF ()
956 : {
957 : gass_type::project_LSF ();
958 : }
959 :
960 : /// checks whether the element is intersected by the interface, or what, and prepares the data
961 : virtual int check_elem_lsf
962 : (
963 : size_t n_co, ///< number of the corners of the element
964 : GridObject * pElem, ///< the element to process
965 : int si, ///< subset of the element
966 : int g_level, ///< grid level of the element
967 : bool use_hanging, ///< if there can be hanging nodes
968 : const MathVector<dim> vCornerCoords [], ///< coordinates of the corners of the element
969 : number time ///< the phisical time
970 : )
971 : {
972 : return gass_type::check_elem_lsf (n_co, pElem, si, g_level, use_hanging, vCornerCoords, time);
973 : }
974 :
975 : /// (slow version) checks whether the element is intersected by the interface, or what, and prepares the data
976 : virtual int check_elem_lsf
977 : (
978 : size_t n_co, ///< number of the corners of the element
979 : GridObject * pElem, ///< the element to process
980 : int si, ///< subset of the element
981 : bool use_hanging, ///< if there can be hanging nodes
982 : const MathVector<dim> vCornerCoords [], ///< coordinates of the corners of the element
983 : number time ///< the phisical time
984 : )
985 : {
986 : ConstSmartPtr<grid_type> mg = base_type::m_spApproxSpace->domain()->grid ();
987 : int g_level = mg->get_level (pElem);
988 : UG_ASSERT (g_level >= 0, "LSGFDomainDiscretization: Grid element without grid level.");
989 : return gass_type::check_elem_lsf (n_co, pElem, si, g_level, use_hanging, vCornerCoords, time);
990 : }
991 :
992 : /// extrapolates a component of the solution to the vertices behind the interface (w.r.t. a base corner)
993 : virtual void extrapolate_by_lsf
994 : (
995 : size_t num_co, ///< number of the corners
996 : size_t base_co, ///< the base corner
997 : number * u, ///< nodal values to extrapolate
998 : size_t fct ///< index of the function (to identify to type of the extrapolation)
999 : ) const
1000 : {
1001 : gass_type::extrapolate_by_lsf (num_co, base_co, u, fct);
1002 : }
1003 :
1004 : /// extrapolates a component of the solution to the vertices behind the interface (by averaging)
1005 : virtual void extrapolate_by_lsf
1006 : (
1007 : size_t num_co, ///< number of the corners
1008 : number * u, ///< nodal values to extrapolate
1009 : size_t fct ///< index of the function (to identify to type of the extrapolation)
1010 : ) const
1011 : {
1012 : gass_type::extrapolate_by_lsf (num_co, u, fct);
1013 : }
1014 :
1015 : /// returns true if the corner is "inside" (use after check_elem_lsf)
1016 : virtual bool corner_inside
1017 : (
1018 : size_t co ///< the corner
1019 : ) const
1020 : {return gass_type::corner_inside (co);}
1021 :
1022 : /// returns the effective value of the LSF at a corner (use after check_elem_lsf)
1023 : virtual number lsf_at
1024 : (
1025 : size_t co ///< the corner
1026 : ) const
1027 : {return gass_type::lsf_at (co);}
1028 :
1029 : protected:
1030 : /// the Level-Set Function constraint
1031 : SmartPtr<ls_constraint_type> m_spLSFGFConstraint;
1032 : };
1033 :
1034 : } // end namespace ug
1035 :
1036 : #include "dom_disc_embb_impl.h"
1037 :
1038 : #endif // __H__UG__PLUGINS__D3F__EMBASS__
1039 :
1040 : /* End of File */
|
