ug4
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"
44 namespace ug{
45 
47 // InverseLinker
49 
50 template <int dim>
52 InverseLinker(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 
64 template <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 
94 template <int dim>
97 {
98  divide(CreateConstUserData<dim>(dividend, number()), divisor);
99 }
100 
101 template <int dim>
103 divide(SmartPtr<CplUserData<number, dim> > dividend, number divisor)
104 {
105  divide(dividend, CreateConstUserData<dim>(divisor, number()));
106 }
107 
108 template <int dim>
110 divide(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
117 template <int dim>
119 evaluate (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 
140 template <int dim>
141 template <int refDim>
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  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 
173 template <int dim>
174 template <int refDim>
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  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:52
#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