Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Ivo Muha, Andreas Vogel
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 : * andreasvogel scale_add_linker_impl.h used as template
35 : */
36 :
37 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__INVERSE_LINKER_IMPL__
38 : #define __H__UG__LIB_DISC__SPATIAL_DISC__INVERSE_LINKER_IMPL__
39 :
40 : #include "inverse_linker.h"
41 : #include "linker_traits.h"
42 : #include "lib_disc/spatial_disc/user_data/const_user_data.h"
43 : #include "common/util/number_util.h"
44 : namespace ug{
45 :
46 : ////////////////////////////////////////////////////////////////////////////////
47 : // InverseLinker
48 : ////////////////////////////////////////////////////////////////////////////////
49 :
50 : template <int dim>
51 0 : InverseLinker<dim>::
52 0 : InverseLinker(const InverseLinker& linker)
53 : {
54 0 : if(linker.m_vpDividendData.size() != linker.m_vpDivisorData.size())
55 0 : UG_THROW("InverseLinker: number of inputs mismatch.");
56 :
57 0 : for(size_t i = 0; i < linker.m_vpDivisorData.size(); ++i)
58 : {
59 0 : this->divide(linker.m_vpDividendData[i], linker.m_vpDivisorData[i]);
60 : }
61 0 : }
62 :
63 :
64 : template <int dim>
65 0 : void InverseLinker<dim>::
66 : divide(SmartPtr<CplUserData<number, dim> > dividend, SmartPtr<CplUserData<number, dim> > divisor)
67 : {
68 :
69 : // current number of inputs
70 0 : const size_t numInput = base_type::num_input() / 2;
71 :
72 : // resize scaling
73 0 : m_vpDivisorData.resize(numInput+1);
74 0 : m_vpDependData.resize(numInput+1);
75 0 : m_vpDividendData.resize(numInput+1);
76 0 : m_vpDividendDependData.resize(numInput+1);
77 :
78 : // remember userdata
79 : UG_ASSERT(divisor.valid(), "Null Pointer as Input set.");
80 0 : m_vpDivisorData[numInput] = divisor;
81 0 : m_vpDependData[numInput] = divisor.template cast_dynamic<DependentUserData<number, dim> >();
82 :
83 : // remember userdata
84 : UG_ASSERT(dividend.valid(), "Null Pointer as Scale set.");
85 0 : m_vpDividendData[numInput] = dividend;
86 0 : m_vpDividendDependData[numInput] = dividend.template cast_dynamic<DependentUserData<number, dim> >();
87 :
88 : // increase number of inputs by one and set inputs at base class
89 0 : base_type::set_num_input(2*numInput+2);
90 0 : base_type::set_input(2*numInput, divisor, divisor);
91 0 : base_type::set_input(2*numInput+1, dividend, dividend);
92 0 : }
93 :
94 : template <int dim>
95 0 : void InverseLinker<dim>::
96 : divide(number dividend, SmartPtr<CplUserData<number, dim> > divisor)
97 : {
98 0 : divide(CreateConstUserData<dim>(dividend, number()), divisor);
99 0 : }
100 :
101 : template <int dim>
102 0 : void InverseLinker<dim>::
103 : divide(SmartPtr<CplUserData<number, dim> > dividend, number divisor)
104 : {
105 0 : divide(dividend, CreateConstUserData<dim>(divisor, number()));
106 0 : }
107 :
108 : template <int dim>
109 0 : void InverseLinker<dim>::
110 : divide(number dividend, number divisor)
111 : {
112 0 : divide(CreateConstUserData<dim>(dividend, number()),
113 : CreateConstUserData<dim>(divisor, number()));
114 0 : }
115 :
116 : //Scale ist Dividend! UserData ist Divisor
117 : template <int dim>
118 0 : void InverseLinker<dim>::
119 : evaluate (number& value,
120 : const MathVector<dim>& globIP,
121 : number time, int si) const
122 : {
123 : // reset value
124 0 : value = 1.0;
125 :
126 0 : number valDivisor=0;
127 0 : number valDividend=0;
128 :
129 0 : for(size_t c = 0; c < m_vpDivisorData.size(); ++c)
130 : {
131 0 : (*m_vpDivisorData[c])(valDivisor, globIP, time, si);
132 0 : (*m_vpDividendData[c])(valDividend, globIP, time, si);
133 : UG_ASSERT(valDivisor!=0, "DIVISOR IS 0");
134 0 : UG_COND_THROW(valDivisor==0, "DIVISOR IS 0");
135 :
136 0 : value *= valDividend/valDivisor;
137 : }
138 0 : }
139 :
140 : template <int dim>
141 : template <int refDim>
142 0 : void InverseLinker<dim>::
143 : evaluate(number vValue[],
144 : const MathVector<dim> vGlobIP[],
145 : number time, int si,
146 : GridObject* elem,
147 : const MathVector<dim> vCornerCoords[],
148 : const MathVector<refDim> vLocIP[],
149 : const size_t nip,
150 : LocalVector* u,
151 : const MathMatrix<refDim, dim>* vJT) const
152 : {
153 : // reset value
154 0 : for(size_t ip = 0; ip < nip; ++ip)
155 0 : vValue[ip] = 1.0;
156 :
157 0 : std::vector<number> vValData(nip);
158 0 : std::vector<number> vValScale(nip);
159 :
160 : // add contribution of each summand
161 0 : for(size_t c = 0; c < m_vpDivisorData.size(); ++c)
162 : {
163 0 : (*m_vpDivisorData[c])(&vValData[0], vGlobIP, time, si,
164 : elem, vCornerCoords, vLocIP, nip, u, vJT);
165 0 : (*m_vpDividendData[c])(&vValScale[0], vGlobIP, time, si,
166 : elem, vCornerCoords, vLocIP, nip, u, vJT);
167 :
168 0 : for(size_t ip = 0; ip < nip; ++ip)
169 0 : vValue[ip] *= vValScale[ip]/vValData[ip];
170 : }
171 0 : }
172 :
173 : template <int dim>
174 : template <int refDim>
175 0 : void InverseLinker<dim>::
176 : eval_and_deriv(number vValue[],
177 : const MathVector<dim> vGlobIP[],
178 : number time, int si,
179 : GridObject* elem,
180 : const MathVector<dim> vCornerCoords[],
181 : const MathVector<refDim> vLocIP[],
182 : const size_t nip,
183 : LocalVector* u,
184 : bool bDeriv,
185 : int s,
186 : std::vector<std::vector<number> > vvvDeriv[],
187 : const MathMatrix<refDim, dim>* vJT) const
188 : {
189 : // check that size of Scalings and inputs is equal
190 : UG_ASSERT(m_vpDivisorData.size() == m_vpDividendData.size(), "Wrong num Scales.");
191 :
192 : // compute value
193 0 : for(size_t ip = 0; ip < nip; ++ip)
194 : {
195 : // reset value
196 0 : vValue[ip] = 1.0;
197 :
198 : // add contribution of each summand
199 0 : for(size_t c2 = 0; c2 < m_vpDivisorData.size(); ++c2)
200 : {
201 0 : UG_COND_THROW(CloseToZero(divisor_value(c2,s,ip)), "DIVISOR IS 0");
202 0 : vValue[ip] *= dividend_value(c2,s,ip)/divisor_value(c2,s,ip);
203 : }
204 : }
205 :
206 : // compute value
207 : // check if derivative is required
208 0 : if(!bDeriv || this->zero_derivative()) return;
209 :
210 : // check sizes
211 : UG_ASSERT(m_vpDependData.size() == m_vpDividendDependData.size(),
212 : "Wrong num Scales.");
213 :
214 : // clear all derivative values
215 0 : this->set_zero(vvvDeriv, nip);
216 :
217 : // (dividend/divisor)' = (u/v)' = (u'v - uv')/v^2 = u'/v - uv'/v^2
218 : //
219 : // now u = prod_i u_i and v = prod_i v_i
220 : // set nu_i = prod_{j != i} u_j and nv_i = prod_{j != i} v_j
221 : // note that prod_i u_i = u_j nu_j for all j.
222 : // then u'_i = sum_j nu_j u_j'
223 : //
224 : // (u/v)' = ((sum_j nu_j u_j') prod_i v_i - (sum_j nv_j v_j') prod_i u_i ) / prod_i v_i^2
225 : // = ((sum_j nu_j u_j' v_j nv_j - sum_j nv_j v_j' u_j nu_j ) / v_j^2 nv_j^2
226 : // = sum_j ( nu_j / nv_j * (u_j' v_j - u_j v_j') / v_j^2 )
227 :
228 :
229 : // loop all inputs
230 0 : for(size_t c = 0; c < m_vpDivisorData.size(); ++c)
231 : {
232 0 : for(size_t ip = 0; ip < nip; ++ip)
233 : {
234 : double vValue = 1.0; // vValue = nu_j / nv_j = prod_{k ! = j} u_k/v_k
235 0 : for(size_t c2 = 0; c2 < m_vpDivisorData.size(); ++c2)
236 : {
237 0 : if(c == c2) continue;
238 0 : vValue *= dividend_value(c2,s,ip)/divisor_value(c2,s,ip);
239 : }
240 :
241 : // check if Divisor has derivative
242 0 : if(!m_vpDivisorData[c]->zero_derivative())
243 : {
244 :
245 : // loop functions
246 0 : for(size_t fct = 0; fct < divisor_num_fct(c); ++fct)
247 : {
248 : // get common fct id for this function
249 : const size_t commonFct = divisor_common_fct(c, fct);
250 :
251 : // loop dofs
252 0 : for(size_t sh = 0; sh < this->num_sh(fct); ++sh)
253 : {
254 0 : if(dividend_value(c,s,ip)!=0)
255 : {
256 :
257 0 : vvvDeriv[ip][commonFct][sh] +=
258 0 : vValue *
259 0 : (-1.0)*
260 0 : (dividend_value(c,s,ip))/divisor_value(c,s,ip)/divisor_value(c,s,ip)
261 0 : *divisor_deriv(c, s, ip, fct, sh);
262 0 : UG_COND_THROW(IsFiniteAndNotTooBig(vvvDeriv[ip][commonFct][sh])==false, "");
263 : }
264 : else
265 0 : vvvDeriv[ip][commonFct][sh] += 0;
266 : // input_deriv(c, s, ip, fct, sh),
267 : // scale_value(c, s, ip));
268 : }
269 : }
270 : }
271 :
272 : // check if Dividend has derivative
273 0 : if(!m_vpDividendData[c]->zero_derivative())
274 : {
275 :
276 : // loop functions
277 0 : for(size_t fct = 0; fct < dividend_num_fct(c); ++fct)
278 : {
279 : // get common fct id for this function
280 : const size_t commonFct = dividend_common_fct(c, fct);
281 :
282 : // loop dofs
283 0 : for(size_t sh = 0; sh < this->num_sh(fct); ++sh)
284 : {
285 0 : if(dividend_value(c,s,ip)!=0)
286 : {
287 0 : vvvDeriv[ip][commonFct][sh] +=
288 : vValue
289 0 : * dividend_deriv(c, s, ip, fct, sh) / divisor_value(c,s,ip);
290 0 : UG_COND_THROW(IsFiniteAndNotTooBig(vvvDeriv[ip][commonFct][sh])==false, "");
291 : }
292 : else
293 0 : vvvDeriv[ip][commonFct][sh] += 0;
294 :
295 : //linker_traits<TData, TDataScale>::
296 : //mult_add(this->deriv(s, ip, commonFct, sh),
297 : // input_value(c, s, ip),
298 : // scale_deriv(c, s, ip, fct, sh));
299 : }
300 : }
301 : }
302 :
303 : }
304 :
305 :
306 : }
307 :
308 : }
309 :
310 :
311 : } // end namespace ug
312 :
313 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__INVERSE_LINKER_IMPL__ */
|
