Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: 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 : #include "data_export.h"
34 : #include "lib_disc/reference_element/reference_element_util.h"
35 :
36 : namespace ug{
37 :
38 : ////////////////////////////////////////////////////////////////////////////////
39 : // ValueDataExport
40 : ////////////////////////////////////////////////////////////////////////////////
41 :
42 : template <int dim>
43 : template <int refDim>
44 0 : void ValueDataExport<dim>::eval_and_deriv(number vValue[],
45 : const MathVector<dim> vGlobIP[],
46 : number time, int si,
47 : GridObject* elem,
48 : const MathVector<dim> vCornerCoords[],
49 : const MathVector<refDim> vLocIP[],
50 : const size_t nip,
51 : LocalVector* u,
52 : bool bDeriv,
53 : int s,
54 : std::vector<std::vector<number> > vvvDeriv[],
55 : const MathMatrix<refDim, dim>* vJT) const
56 : {
57 : // abbreviation for component
58 : static const int _C_ = 0;
59 :
60 : // reference object id
61 0 : const ReferenceObjectID roid = elem->reference_object_id();
62 :
63 : // local finite element id
64 0 : const LFEID& lfeID = this->function_group().local_finite_element_id(_C_);
65 :
66 : // access local vector by map
67 0 : u->access_by_map(this->map());
68 :
69 : // request for trial space
70 : try{
71 : const LocalShapeFunctionSet<refDim>& rTrialSpace
72 : = LocalFiniteElementProvider::get<refDim>(roid, lfeID);
73 :
74 : // memory for shapes
75 : std::vector<number> vShape;
76 :
77 : // loop ips
78 0 : for(size_t ip = 0; ip < nip; ++ip)
79 : {
80 : // evaluate at shapes at ip
81 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
82 :
83 : // compute value at ip
84 0 : vValue[ip] = 0.0;
85 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
86 0 : vValue[ip] += (*u)(_C_, sh) * vShape[sh];
87 :
88 : // store derivative
89 0 : if(bDeriv)
90 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
91 0 : vvvDeriv[ip][_C_][sh] = vShape[sh];
92 : }
93 :
94 0 : }
95 0 : UG_CATCH_THROW("ValueDataExport: Trial space missing, Reference Object: "
96 : <<roid<<", Trial Space: "<<lfeID<<", refDim="<<refDim);
97 0 : }
98 :
99 : template <int dim>
100 0 : void ValueDataExport<dim>::check_setup() const
101 : {
102 0 : if((int)this->function_group().size() != 1)
103 0 : UG_THROW("ValueDataExport: Expected to work on "<<1<<" function, "
104 : "but only set "<<this->function_group().size()<<" ("<<
105 : this->function_group().names()<<").");
106 0 : }
107 :
108 : template <int dim>
109 0 : bool ValueDataExport<dim>::continuous() const
110 : {
111 : static const int _C_ = 0;
112 0 : const LFEID& lfeID = this->function_group().local_finite_element_id(_C_);
113 0 : return LocalFiniteElementProvider::continuous(lfeID);
114 : }
115 :
116 : ////////////////////////////////////////////////////////////////////////////////
117 : // GradientDataExport
118 : ////////////////////////////////////////////////////////////////////////////////
119 :
120 : template <int dim>
121 : template <int refDim>
122 0 : void GradientDataExport<dim>::eval_and_deriv(MathVector<dim> vValue[],
123 : const MathVector<dim> vGlobIP[],
124 : number time, int si,
125 : GridObject* elem,
126 : const MathVector<dim> vCornerCoords[],
127 : const MathVector<refDim> vLocIP[],
128 : const size_t nip,
129 : LocalVector* u,
130 : bool bDeriv,
131 : int s,
132 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[],
133 : const MathMatrix<refDim, dim>* vJT) const
134 : {
135 : // abbreviation for component
136 : static const int _C_ = 0;
137 :
138 : // reference object id
139 0 : const ReferenceObjectID roid = elem->reference_object_id();
140 :
141 : // local finite element id
142 0 : const LFEID& lfeID = this->function_group().local_finite_element_id(_C_);
143 :
144 : // access local vector by map
145 0 : u->access_by_map(this->map());
146 :
147 : // request for trial space
148 : try{
149 : const LocalShapeFunctionSet<refDim>& rTrialSpace
150 : = LocalFiniteElementProvider::get<refDim>(roid, lfeID);
151 :
152 : // Reference Mapping
153 : MathMatrix<dim, refDim> JTInv;
154 : std::vector<MathMatrix<refDim, dim> > vJTtmp;
155 0 : if(!vJT){
156 : DimReferenceMapping<refDim, dim>& map
157 : = ReferenceMappingProvider::get<refDim, dim>(roid, vCornerCoords);
158 :
159 0 : vJTtmp.resize(nip);
160 0 : map.jacobian_transposed(&vJTtmp[0], vLocIP, nip);
161 : vJT = &vJTtmp[0];
162 : }
163 :
164 : // storage for shape function at ip
165 : std::vector<MathVector<refDim> > vLocGrad;
166 : MathVector<refDim> locGrad;
167 :
168 : // loop ips
169 0 : for(size_t ip = 0; ip < nip; ++ip)
170 : {
171 : // evaluate at shapes at ip
172 0 : rTrialSpace.grads(vLocGrad, vLocIP[ip]);
173 :
174 : // compute grad at ip
175 : VecSet(locGrad, 0.0);
176 0 : for(size_t sh = 0; sh < vLocGrad.size(); ++sh)
177 0 : VecScaleAppend(locGrad, (*u)(_C_, sh), vLocGrad[sh]);
178 :
179 0 : Inverse(JTInv, vJT[ip]);
180 0 : MatVecMult(vValue[ip], JTInv, locGrad);
181 :
182 : // store derivative
183 0 : if(bDeriv)
184 0 : for(size_t sh = 0; sh < vLocGrad.size(); ++sh)
185 0 : MatVecMult(vvvDeriv[ip][_C_][sh], JTInv, vLocGrad[sh]);
186 : }
187 :
188 0 : }
189 0 : UG_CATCH_THROW("GradientDataExport: Trial space missing, Reference Object: "
190 : <<roid<<", Trial Space: "<<lfeID<<", refDim="<<refDim);
191 0 : }
192 :
193 : template <int dim>
194 0 : void GradientDataExport<dim>::check_setup() const
195 : {
196 0 : if((int)this->function_group().size() != 1)
197 0 : UG_THROW("GradientDataExport: Expected to work on "<<1<<" function, "
198 : "but only set "<<this->function_group().size()<<" ("<<
199 : this->function_group().names()<<").");
200 0 : }
201 :
202 : ////////////////////////////////////////////////////////////////////////////////
203 : // VectorDataExport
204 : ////////////////////////////////////////////////////////////////////////////////
205 :
206 : template <int dim>
207 : template <int refDim>
208 0 : void VectorDataExport<dim>::eval_and_deriv(MathVector<dim> vValue[],
209 : const MathVector<dim> vGlobIP[],
210 : number time, int si,
211 : GridObject* elem,
212 : const MathVector<dim> vCornerCoords[],
213 : const MathVector<refDim> vLocIP[],
214 : const size_t nip,
215 : LocalVector* u,
216 : bool bDeriv,
217 : int s,
218 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[],
219 : const MathMatrix<refDim, dim>* vJT) const
220 : {
221 : // reference object id
222 0 : const ReferenceObjectID roid = elem->reference_object_id();
223 :
224 : // access local vector by map
225 0 : u->access_by_map(this->map());
226 :
227 0 : if(bDeriv)
228 0 : this->set_zero(vvvDeriv, nip);
229 :
230 0 : for(size_t d = 0; d < dim; ++d)
231 : {
232 : // local finite element id
233 0 : const LFEID& lfeID = this->function_group().local_finite_element_id(d);
234 :
235 : // request for trial space
236 : try{
237 : const LocalShapeFunctionSet<refDim>& rTrialSpace
238 : = LocalFiniteElementProvider::get<refDim>(roid, lfeID);
239 :
240 : // memory for shapes
241 : std::vector<number> vShape;
242 :
243 : // loop ips
244 0 : for(size_t ip = 0; ip < nip; ++ip)
245 : {
246 : // evaluate at shapes at ip
247 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
248 :
249 : // compute value at ip
250 0 : vValue[ip] = 0.0;
251 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
252 0 : vValue[ip][d] += (*u)(d, sh) * vShape[sh];
253 :
254 : // store derivative
255 0 : if(bDeriv)
256 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
257 0 : vvvDeriv[ip][d][sh][d] = vShape[sh];
258 : }
259 :
260 0 : }
261 0 : UG_CATCH_THROW("VectorDataExport: Trial space missing, Reference Object: "
262 : <<roid<<", Trial Space: "<<lfeID<<", refDim="<<refDim);
263 : }
264 0 : }
265 :
266 : template <int dim>
267 0 : void VectorDataExport<dim>::check_setup() const
268 : {
269 0 : if((int)this->function_group().size() != dim)
270 0 : UG_THROW("VectorDataExport: Expected to work on "<<dim<<" functions, "
271 : "but only set "<<this->function_group().size()<<" ("<<
272 : this->function_group().names()<<").");
273 0 : }
274 :
275 : template <int dim>
276 0 : bool VectorDataExport<dim>::continuous() const
277 : {
278 0 : for(size_t d = 0; d < dim; ++d){
279 0 : const LFEID& lfeID = this->function_group().local_finite_element_id(d);
280 0 : if(!LocalFiniteElementProvider::continuous(lfeID)) return false;
281 : }
282 : return true;
283 : }
284 :
285 : ////////////////////////////////////////////////////////////////////////////////
286 : // explicit instantiations
287 : ////////////////////////////////////////////////////////////////////////////////
288 :
289 : #ifdef UG_DIM_1
290 : template class ValueDataExport<1>;
291 : template class GradientDataExport<1>;
292 : template class VectorDataExport<1>;
293 : #endif
294 : #ifdef UG_DIM_2
295 : template class ValueDataExport<2>;
296 : template class GradientDataExport<2>;
297 : template class VectorDataExport<2>;
298 : #endif
299 : #ifdef UG_DIM_3
300 : template class ValueDataExport<3>;
301 : template class GradientDataExport<3>;
302 : template class VectorDataExport<3>;
303 : #endif
304 :
305 : } // end namespace ug
|
