Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
inverse_linker_impl.h
Go to the documentation of this file.
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"
44namespace ug{
45
47// InverseLinker
49
50template <int dim>
52InverseLinker(const InverseLinker& linker)
53{
54 if(linker.m_vpDividendData.size() != linker.m_vpDivisorData.size())
55 UG_THROW("InverseLinker: number of inputs mismatch.");
56
57 for(size_t i = 0; i < linker.m_vpDivisorData.size(); ++i)
58 {
59 this->divide(linker.m_vpDividendData[i], linker.m_vpDivisorData[i]);
60 }
61}
62
63
64template <int dim>
67{
68
69// current number of inputs
70 const size_t numInput = base_type::num_input() / 2;
71
72// resize scaling
73 m_vpDivisorData.resize(numInput+1);
74 m_vpDependData.resize(numInput+1);
75 m_vpDividendData.resize(numInput+1);
76 m_vpDividendDependData.resize(numInput+1);
77
78// remember userdata
79 UG_ASSERT(divisor.valid(), "Null Pointer as Input set.");
80 m_vpDivisorData[numInput] = divisor;
81 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 m_vpDividendData[numInput] = dividend;
86 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 base_type::set_num_input(2*numInput+2);
90 base_type::set_input(2*numInput, divisor, divisor);
91 base_type::set_input(2*numInput+1, dividend, dividend);
92}
93
94template <int dim>
97{
98 divide(CreateConstUserData<dim>(dividend, number()), divisor);
99}
100
101template <int dim>
104{
105 divide(dividend, CreateConstUserData<dim>(divisor, number()));
106}
107
108template <int dim>
110divide(number dividend, number divisor)
111{
112 divide(CreateConstUserData<dim>(dividend, number()),
113 CreateConstUserData<dim>(divisor, number()));
114}
115
116//Scale ist Dividend! UserData ist Divisor
117template <int dim>
119evaluate (number& value,
120 const MathVector<dim>& globIP,
121 number time, int si) const
122{
123 // reset value
124 value = 1.0;
125
126 number valDivisor=0;
127 number valDividend=0;
128
129 for(size_t c = 0; c < m_vpDivisorData.size(); ++c)
130 {
131 (*m_vpDivisorData[c])(valDivisor, globIP, time, si);
132 (*m_vpDividendData[c])(valDividend, globIP, time, si);
133 UG_ASSERT(valDivisor!=0, "DIVISOR IS 0");
134 UG_COND_THROW(valDivisor==0, "DIVISOR IS 0");
135
136 value *= valDividend/valDivisor;
137 }
138}
139
140template <int dim>
141template <int refDim>
143evaluate(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 for(size_t ip = 0; ip < nip; ++ip)
155 vValue[ip] = 1.0;
156
157 std::vector<number> vValData(nip);
158 std::vector<number> vValScale(nip);
159
160// add contribution of each summand
161 for(size_t c = 0; c < m_vpDivisorData.size(); ++c)
162 {
163 (*m_vpDivisorData[c])(&vValData[0], vGlobIP, time, si,
164 elem, vCornerCoords, vLocIP, nip, u, vJT);
165 (*m_vpDividendData[c])(&vValScale[0], vGlobIP, time, si,
166 elem, vCornerCoords, vLocIP, nip, u, vJT);
167
168 for(size_t ip = 0; ip < nip; ++ip)
169 vValue[ip] *= vValScale[ip]/vValData[ip];
170 }
171}
172
173template <int dim>
174template <int refDim>
176eval_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 for(size_t ip = 0; ip < nip; ++ip)
194 {
195 // reset value
196 vValue[ip] = 1.0;
197
198 // add contribution of each summand
199 for(size_t c2 = 0; c2 < m_vpDivisorData.size(); ++c2)
200 {
201 UG_COND_THROW(CloseToZero(divisor_value(c2,s,ip)), "DIVISOR IS 0");
202 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 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 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 for(size_t c = 0; c < m_vpDivisorData.size(); ++c)
231 {
232 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 for(size_t c2 = 0; c2 < m_vpDivisorData.size(); ++c2)
236 {
237 if(c == c2) continue;
238 vValue *= dividend_value(c2,s,ip)/divisor_value(c2,s,ip);
239 }
240
241 // check if Divisor has derivative
242 if(!m_vpDivisorData[c]->zero_derivative())
243 {
244
245 // loop functions
246 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 for(size_t sh = 0; sh < this->num_sh(fct); ++sh)
253 {
254 if(dividend_value(c,s,ip)!=0)
255 {
256
257 vvvDeriv[ip][commonFct][sh] +=
258 vValue *
259 (-1.0)*
260 (dividend_value(c,s,ip))/divisor_value(c,s,ip)/divisor_value(c,s,ip)
261 *divisor_deriv(c, s, ip, fct, sh);
262 UG_COND_THROW(IsFiniteAndNotTooBig(vvvDeriv[ip][commonFct][sh])==false, "");
263 }
264 else
265 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 if(!m_vpDividendData[c]->zero_derivative())
274 {
275
276 // loop functions
277 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 for(size_t sh = 0; sh < this->num_sh(fct); ++sh)
284 {
285 if(dividend_value(c,s,ip)!=0)
286 {
287 vvvDeriv[ip][commonFct][sh] +=
288 vValue
289 * dividend_deriv(c, s, ip, fct, sh) / divisor_value(c,s,ip);
290 UG_COND_THROW(IsFiniteAndNotTooBig(vvvDeriv[ip][commonFct][sh])==false, "");
291 }
292 else
293 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__ */
parameterString s
Definition smart_pointer.h:108
Type based UserData.
Definition user_data.h:501
The base class for all geometric objects, such as vertices, edges, faces, volumes,...
Definition grid_base_objects.h:157
Definition inverse_linker.h:60
InverseLinker()
constructor
Definition inverse_linker.h:67
void evaluate(number &value, const MathVector< dim > &globIP, number time, int si) const
Definition inverse_linker_impl.h:119
std::vector< SmartPtr< CplUserData< number, dim > > > m_vpDivisorData
data Divisor
Definition inverse_linker.h:163
void eval_and_deriv(number vValue[], const MathVector< dim > vGlobIP[], number time, int si, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< refDim > vLocIP[], const size_t nip, LocalVector *u, bool bDeriv, int s, std::vector< std::vector< number > > vvvDeriv[], const MathMatrix< refDim, dim > *vJT=NULL) const
Definition inverse_linker_impl.h:176
void divide(SmartPtr< CplUserData< number, dim > > dividend, SmartPtr< CplUserData< number, dim > > divisor)
Definition inverse_linker_impl.h:66
std::vector< SmartPtr< CplUserData< number, dim > > > m_vpDividendData
data Dividend
Definition inverse_linker.h:157
Definition local_algebra.h:198
A class for fixed size, dense matrices.
Definition math_matrix.h:63
a mathematical Vector with N entries.
Definition math_vector.h:97
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define UG_THROW(msg)
Definition error.h:57
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
double number
Definition types.h:124
the ug namespace
bool CloseToZero(double d)
Definition number_util.h:48
bool IsFiniteAndNotTooBig(double d)
Definition number_util.h:39