Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
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 : #ifndef __H__LIB_ALGEBRA__OPERATOR__INTERFACE__CONSTRAINED_LINEAR_ITERATOR__
34 : #define __H__LIB_ALGEBRA__OPERATOR__INTERFACE__CONSTRAINED_LINEAR_ITERATOR__
35 :
36 : #include "common/util/smart_pointer.h"
37 : #include "lib_algebra/operator/interface/linear_iterator.h"
38 : #include "lib_disc/function_spaces/grid_function.h"
39 : #include "lib_disc/spatial_disc/domain_disc_interface.h"
40 : #include "preconditioner.h"
41 :
42 : #include <type_traits>
43 :
44 : namespace ug {
45 :
46 :
47 : /**
48 : * This class is a template for constraint-respecting versions of
49 : * any linear iterator derived from ILinearIterator. This includes
50 : * linear solvers as well as preconditioners.
51 : * The derived ILinearIterator class must have a constructor
52 : * without arguments.
53 : * The handled constraints must implement the two methods:
54 : * adjust_correction(),
55 : * adjust_defect().
56 : *
57 : * @note This class may have lost its use:
58 : * Constraints are now usually implemented in such a way that
59 : * that both defect and correction is always zero in constrained
60 : * DoFs.
61 : */
62 : template <typename TDomain, typename TAlgebra, typename TLinIt, typename = void>
63 : class ConstrainedLinearIterator : public TLinIt
64 : {
65 : public:
66 : typedef GridFunction<TDomain, TAlgebra> gf_type;
67 : typedef typename TAlgebra::vector_type vector_type;
68 : typedef TLinIt base_type;
69 :
70 : using TLinIt::init;
71 :
72 : public:
73 : /// constructor
74 0 : ConstrainedLinearIterator(SmartPtr<IDomainDiscretization<TAlgebra> > domDisc)
75 0 : : base_type(), m_spDomDisc(domDisc), m_time(0.0)
76 0 : {}
77 :
78 : /// clone constructor
79 0 : ConstrainedLinearIterator(const ConstrainedLinearIterator<TDomain, TAlgebra, TLinIt> &parent)
80 0 : : base_type(parent), m_spDomDisc(parent.m_spDomDisc), m_time(parent.m_time)
81 0 : {}
82 :
83 : /// clone
84 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
85 : {
86 0 : return make_sp(new ConstrainedLinearIterator<TDomain, TAlgebra, TLinIt>(*this));
87 : }
88 :
89 : /// destructor
90 0 : virtual ~ConstrainedLinearIterator(){}
91 :
92 : protected:
93 : /// @copydoc ILinearIterator::name()
94 0 : virtual const char* name() const
95 : {
96 0 : return base_type::name();
97 : }
98 :
99 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type,vector_type> > J, const vector_type& u)
100 : {
101 0 : return base_type::init(J, u);
102 : }
103 :
104 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type,vector_type> > L)
105 : {
106 0 : return base_type::init(L);
107 : }
108 :
109 0 : virtual bool apply(vector_type& c, const vector_type& d)
110 : {
111 : // perform normal step
112 0 : if (!base_type::apply(c, d)) return false;
113 :
114 : // apply constraints on correction
115 0 : if (!m_spDomDisc.valid())
116 0 : UG_THROW("Domain discretization passed to ConstrainedLinearIterator is not valid.\n"
117 : "Cannot apply any constraints.")
118 :
119 0 : gf_type* gf = dynamic_cast<gf_type*>(&c);
120 0 : if (gf && gf->dof_distribution().valid())
121 : {
122 0 : size_t nConstr = m_spDomDisc->num_constraints();
123 :
124 : // first Dirichlet (hanging nodes might be constrained by Dirichlet nodes)
125 0 : for (size_t i = 0; i < nConstr; i++)
126 0 : if (m_spDomDisc->constraint(i)->type() & CT_DIRICHLET)
127 0 : m_spDomDisc->constraint(i)->adjust_correction(c, gf->dof_distribution(), CT_DIRICHLET, m_time);
128 :
129 : // then other interpolation
130 0 : for (size_t i = 0; i < nConstr; i++)
131 0 : if (m_spDomDisc->constraint(i)->type() & CT_CONSTRAINTS)
132 0 : m_spDomDisc->constraint(i)->adjust_correction(c, gf->dof_distribution(), CT_CONSTRAINTS, m_time);
133 :
134 : // then hanging nodes
135 0 : for (size_t i = 0; i < nConstr; i++)
136 0 : if (m_spDomDisc->constraint(i)->type() & CT_HANGING)
137 0 : m_spDomDisc->constraint(i)->adjust_correction(c, gf->dof_distribution(), CT_HANGING, m_time);
138 :
139 : // and Dirichlet again (Dirichlet nodes might also be hanging)
140 0 : for (size_t i = 0; i < nConstr; i++)
141 0 : if (m_spDomDisc->constraint(i)->type() & CT_DIRICHLET)
142 0 : m_spDomDisc->constraint(i)->adjust_correction(c, gf->dof_distribution(), CT_DIRICHLET, m_time);
143 : }
144 0 : else UG_THROW("Vector is not a grid function. This is not supported.")
145 :
146 : return true;
147 : }
148 :
149 0 : virtual bool apply_update_defect(vector_type& c, vector_type& d)
150 : {
151 : apply_update_defect_impl<TLinIt> impl(*this);
152 0 : return impl(c, d);
153 : }
154 :
155 : private:
156 : template <typename S, typename = void>//, typename otherDummy = void>
157 : struct apply_update_defect_impl
158 : {
159 : ConstrainedLinearIterator& cli;
160 :
161 0 : explicit apply_update_defect_impl(ConstrainedLinearIterator& _cli)
162 0 : : cli(_cli) {}
163 :
164 0 : bool operator()(vector_type& c, vector_type& d)
165 : {
166 : // FIXME: This implementation is not correct.
167 : // As the defect is calculated from possibly erroneous corrections, there is no telling
168 : // how to adjust the defect. Except for setting it to zero everywhere, which, by the way,
169 : // is NOT what adjust_defect does, in general. We would need something like
170 : // adjust_defect_linearSolver() here to distinguish between the non-linear defect and
171 : // the residual in the linear solver.
172 : // We CANNOT use our own apply() here and then update the defect ourselves as the
173 : // ILinearIterator does not have an update method nor any way of applying the underlying
174 : // linear operator.
175 :
176 : // perform normal step
177 0 : if (!cli.base_type::apply_update_defect(c, d)) return false;
178 :
179 : // apply constraints on correction and defect
180 0 : if (!cli.m_spDomDisc.valid())
181 0 : UG_THROW("Domain discretization passed to ConstrainedLinearIterator is not valid.\n"
182 : "Cannot apply any constraints.")
183 :
184 0 : gf_type* gf = dynamic_cast<gf_type*>(&c);
185 0 : if (gf && gf->dof_distribution().valid())
186 : {
187 0 : size_t nConstr = cli.m_spDomDisc->num_constraints();
188 :
189 : // corrections
190 : // first Dirichlet (hanging nodes might be constrained by Dirichlet nodes)
191 0 : for (size_t i = 0; i < nConstr; i++)
192 0 : if (cli.m_spDomDisc->constraint(i)->type() & CT_DIRICHLET)
193 0 : cli.m_spDomDisc->constraint(i)->adjust_correction(c, gf->dof_distribution(), CT_DIRICHLET, cli.m_time);
194 :
195 : // then other interpolation
196 0 : for (size_t i = 0; i < nConstr; i++)
197 0 : if (cli.m_spDomDisc->constraint(i)->type() & CT_CONSTRAINTS)
198 0 : cli.m_spDomDisc->constraint(i)->adjust_correction(c, gf->dof_distribution(), CT_CONSTRAINTS, cli.m_time);
199 :
200 : // then hanging nodes
201 0 : for (size_t i = 0; i < nConstr; i++)
202 0 : if (cli.m_spDomDisc->constraint(i)->type() & CT_HANGING)
203 0 : cli.m_spDomDisc->constraint(i)->adjust_correction(c, gf->dof_distribution(), CT_HANGING, cli.m_time);
204 :
205 : // and Dirichlet again (Dirichlet nodes might also be hanging)
206 0 : for (size_t i = 0; i < nConstr; i++)
207 0 : if (cli.m_spDomDisc->constraint(i)->type() & CT_DIRICHLET)
208 0 : cli.m_spDomDisc->constraint(i)->adjust_correction(c, gf->dof_distribution(), CT_DIRICHLET, cli.m_time);
209 :
210 : // residuals
211 : // hanging nodes first
212 0 : for (size_t i = 0; i < nConstr; i++)
213 0 : if (cli.m_spDomDisc->constraint(i)->type() & CT_HANGING)
214 0 : cli.m_spDomDisc->constraint(i)->adjust_linear_residual(d, c, gf->dof_distribution(), CT_HANGING, cli.m_time);
215 :
216 : // then other constraints
217 0 : for (size_t i = 0; i < nConstr; i++)
218 0 : if (cli.m_spDomDisc->constraint(i)->type() & CT_CONSTRAINTS)
219 0 : cli.m_spDomDisc->constraint(i)->adjust_linear_residual(d, c, gf->dof_distribution(), CT_CONSTRAINTS, cli.m_time);
220 :
221 : // and Dirichlet last
222 0 : for (size_t i = 0; i < nConstr; i++)
223 0 : if (cli.m_spDomDisc->constraint(i)->type() & CT_DIRICHLET)
224 0 : cli.m_spDomDisc->constraint(i)->adjust_linear_residual(d, c, gf->dof_distribution(), CT_DIRICHLET, cli.m_time);
225 :
226 : }
227 0 : else UG_THROW("Vector is not a grid function. This is not supported.")
228 :
229 : return true;
230 : }
231 : };
232 :
233 :
234 : // special implementation for IPreconditioner
235 : template <typename S>
236 : struct apply_update_defect_impl<S,
237 : typename std::enable_if<std::is_base_of<IPreconditioner<TAlgebra>, S>::value>::type>
238 : {
239 : ConstrainedLinearIterator& cli;
240 :
241 : explicit apply_update_defect_impl(ConstrainedLinearIterator& _cli)
242 : : cli(_cli) {}
243 :
244 : bool operator()(vector_type& c, vector_type& d)
245 : {
246 : // apply precond and constraints
247 : if (!cli.apply(c, d)) return false;
248 :
249 : // calculate new defect from preconditioner defect operator
250 : cli.m_spDefectOperator->apply_sub(d, c);
251 :
252 : return true;
253 : }
254 : };
255 :
256 : friend struct apply_update_defect_impl<TLinIt>;
257 :
258 :
259 : public:
260 : /// @copydoc ILinearIterator::supports_parallel()
261 0 : virtual bool supports_parallel() const
262 : {
263 0 : return base_type::supports_parallel();
264 : }
265 :
266 : /// setter for time
267 0 : void set_time(number time)
268 : {
269 0 : m_time = time;
270 0 : };
271 :
272 : protected:
273 : SmartPtr<ILinearIterator<vector_type> > m_spLinIt;
274 : SmartPtr<IDomainDiscretization<TAlgebra> > m_spDomDisc;
275 : number m_time;
276 : };
277 :
278 :
279 :
280 :
281 : } // end namespace ug
282 :
283 :
284 : #endif // __H__LIB_ALGEBRA__OPERATOR__INTERFACE__CONSTRAINED_LINEAR_ITERATOR__
|
