Line data Source code
1 : /*
2 : * Copyright (c) 2011-2017: 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 <sstream>
34 :
35 : #include "data_evaluator.h"
36 : #include "lib_disc/common/groups_util.h"
37 :
38 : namespace ug{
39 :
40 : DebugID DID_DATA_EVALUATOR("DATA_EVALUATOR");
41 :
42 : ///////////////////////////////////////////////////////////////////////////////
43 : // prepare / finish
44 : ///////////////////////////////////////////////////////////////////////////////
45 :
46 : template <typename TDomain>
47 0 : void DataEvaluator<TDomain>::
48 : prepare_elem_loop(const ReferenceObjectID id, int si)
49 : {
50 : // prepare loop (elem disc set local ip series here)
51 : try{
52 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
53 0 : m_vElemDisc[PT_ALL][i]->do_prep_elem_loop(id, si);
54 : }
55 0 : UG_CATCH_THROW("DataEvaluatorBase::prepare_elem_loop: "
56 : "Cannot prepare element loop.");
57 :
58 : // extract data imports and userdatas
59 : try{
60 0 : extract_imports_and_userdata(si, m_discPart);
61 : }
62 0 : UG_CATCH_THROW("DataEvaluatorBase::prepare_elem_loop: "
63 : "Cannot extract imports and userdata.");
64 :
65 : // check setup of imports
66 : try{
67 0 : for(size_t i = 0; i < m_vImport[PT_ALL][MASS].size(); ++i)
68 0 : m_vImport[PT_ALL][MASS][i]->check_setup();
69 0 : for(size_t i = 0; i < m_vImport[PT_ALL][STIFF].size(); ++i)
70 0 : m_vImport[PT_ALL][STIFF][i]->check_setup();
71 0 : for(size_t i = 0; i < m_vImport[PT_ALL][RHS].size(); ++i)
72 0 : m_vImport[PT_ALL][RHS][i]->check_setup();
73 : }
74 0 : UG_CATCH_THROW("DataEvaluatorBase::prepare_elem_loop: Import correctly implemented.");
75 :
76 : // prepare and check dependent data
77 : try{
78 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i){
79 0 : m_vDependentData[i]->check_setup();
80 : }
81 : }
82 0 : UG_CATCH_THROW("DataEvaluatorBase::prepare_elem_loop: Dependent UserData "
83 : " (e.g. Linker or Export) is not ready for evaluation.");
84 :
85 : // evaluate constant data
86 0 : for(size_t i = 0; i < m_vConstData.size(); ++i)
87 0 : m_vConstData[i]->compute((LocalVector*)NULL, NULL, NULL, false);
88 0 : }
89 :
90 : template <typename TDomain>
91 0 : void DataEvaluator<TDomain>::finish_elem_loop()
92 : {
93 : // finish each elem disc
94 : try{
95 0 : for(size_t d = 0; d < m_vElemDisc[PT_ALL].size(); ++d)
96 : {
97 0 : IElemDisc<TDomain>* disc = m_vElemDisc[PT_ALL][d];
98 :
99 0 : disc->do_fsh_elem_loop();
100 :
101 : /* TODO:
102 : * In prepare_elem_loop, the elemdiscs initialize the local ip's independently
103 : * on if they are used. For ex., the ip's used only for the mass matrix are
104 : * initialized, too, even if only the stiffness part is assembled. These ip's
105 : * are not cleared below as they do not get into the lists, and this creates
106 : * issues with the linkers that share subordinated userdata objects. For that,
107 : * we clear here all the assigned ip's.
108 : *
109 : * Should it be done here on in do_fsh_elem_loop?
110 : */
111 0 : for (size_t i = 0; i < disc->num_imports(); ++i)
112 : {
113 : IDataImport<dim>& imp = disc->get_import(i);
114 0 : if(imp.data_given())
115 0 : imp.data()->clear();
116 : }
117 : }
118 : }
119 0 : UG_CATCH_THROW("DataEvaluatorBase::fsh_elem_loop: Cannot finish element loop");
120 :
121 : // clear positions at user data
122 : /* TODO:
123 : * Could it be done in a more elegant way? For ex., why clearing here all the ip
124 : * series and not only the ones assigned to the particular userdata objects?
125 : */
126 0 : clear_positions_in_user_data();
127 0 : }
128 :
129 :
130 : ///////////////////////////////////////////////////////////////////////////////
131 : // Assemble routines
132 : ///////////////////////////////////////////////////////////////////////////////
133 :
134 : template <typename TDomain>
135 0 : void DataEvaluator<TDomain>::prepare_timestep(number future_time, const number time, VectorProxyBase* u, size_t algebra_id)
136 : {
137 : try
138 : {
139 0 : for (size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
140 0 : m_vElemDisc[PT_ALL][i]->do_prep_timestep(future_time, time, u, algebra_id);
141 : }
142 0 : UG_CATCH_THROW("DataEvaluatorBase::prep_timestep: Cannot prepare time step.");
143 0 : }
144 :
145 : template <typename TDomain>
146 0 : void DataEvaluator<TDomain>::
147 : prepare_timestep_elem(const number time, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
148 : {
149 : try{
150 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
151 0 : m_vElemDisc[PT_ALL][i]->do_prep_timestep_elem(time, u, elem, vCornerCoords);
152 : }
153 0 : UG_CATCH_THROW("DataEvaluatorBase::prep_timestep_elem: Cannot prepare timestep.");
154 0 : }
155 :
156 :
157 : template <typename TDomain>
158 0 : void DataEvaluator<TDomain>::
159 : prepare_elem(LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[],
160 : const LocalIndices& ind,
161 : bool bDeriv)
162 : {
163 : // prepare element
164 : try
165 : {
166 : // NOTE: constant data is not processed, since constant == independent of the element
167 0 : for(size_t i = 0; i < m_vPosData.size(); ++i)
168 0 : m_vPosData[i]->prepare_element(elem, vCornerCoords);
169 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i)
170 0 : m_vDependentData[i]->prepare_element(elem, vCornerCoords);
171 :
172 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i){
173 : UG_DLOG(DID_DATA_EVALUATOR, 2, ">>OCT_DISC_DEBUG: " << "data_evaluator.cpp: " << "DataEvaluatorBase.prepare_elem(): m_vElemDisc[PT_ALL][i]->do_prep_elem() " << roid << std::endl);
174 0 : m_vElemDisc[PT_ALL][i]->do_prep_elem(u, elem, roid, vCornerCoords);
175 : }
176 : }
177 0 : UG_CATCH_THROW("DataEvaluatorBase::prep_elem: Cannot prepare element.");
178 :
179 : // adjust lin defect array of imports and derivative array of exports
180 : // INFO: This is place here, since the 'prepare_elem' method of an element
181 : // disc may change the number of integration points, even if the type
182 : // of the element (e.g. triangle, quad) stays the same. This is the
183 : // case for, e.g., the NeumannBoundary element disc.
184 0 : if(bDeriv)
185 : {
186 0 : for(size_t i = 0; i < m_vImport[PT_ALL][MASS].size(); ++i)
187 0 : m_vImport[PT_ALL][MASS][i]->update_dof_sizes(ind);
188 0 : for(size_t i = 0; i < m_vImport[PT_ALL][STIFF].size(); ++i)
189 0 : m_vImport[PT_ALL][STIFF][i]->update_dof_sizes(ind);
190 0 : for(size_t i = 0; i < m_vImport[PT_ALL][RHS].size(); ++i)
191 0 : m_vImport[PT_ALL][RHS][i]->update_dof_sizes(ind);
192 :
193 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i)
194 0 : m_vDependentData[i]->update_dof_sizes(ind);
195 : }
196 :
197 : // evaluate position data
198 0 : for(size_t i = 0; i < m_vPosData.size(); ++i)
199 0 : m_vPosData[i]->compute(&u, elem, vCornerCoords, false);
200 :
201 : // process dependent data:
202 : // We can not simply compute exports first, then Linker, because an export
203 : // itself could depend on other data if implemented somehow in the IElemDisc
204 : // (e.g. using data from some DataImport). Thus, we have to loop the sorted
205 : // vector of all dependent data (that is correctly sorted the way that always
206 : // needed data has previously computed).
207 :
208 : // compute the data
209 : try{
210 0 : if (! time_series_needed ()) { // assemble for the given LocalVector
211 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i){
212 0 : u.access_by_map(m_vDependentData[i]->map());
213 0 : m_vDependentData[i]->compute(&u, elem, vCornerCoords, bDeriv);
214 : }
215 : }
216 : else { // assemble for LocalVectorTimeSeries
217 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i){
218 0 : u.access_by_map(m_vDependentData[i]->map());
219 0 : m_vDependentData[i]->compute(m_pLocTimeSeries, elem, vCornerCoords, bDeriv);
220 : }
221 : }
222 : }
223 0 : UG_CATCH_THROW("DataEvaluatorBase::prep_elem: Cannot compute data for Export or Linker.");
224 0 : }
225 :
226 : template <typename TDomain>
227 0 : void DataEvaluator<TDomain>::finish_timestep(const number time, VectorProxyBase* u, size_t algebra_id)
228 : {
229 : try
230 : {
231 0 : for (size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
232 0 : m_vElemDisc[PT_ALL][i]->do_fsh_timestep(time, u, algebra_id);
233 : }
234 0 : UG_CATCH_THROW("DataEvaluatorBase::finish_timestep: Cannot prepare time step.");
235 0 : }
236 :
237 : template <typename TDomain>
238 0 : void DataEvaluator<TDomain>::
239 : finish_timestep_elem(const number time, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
240 : {
241 : try{
242 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
243 0 : m_vElemDisc[PT_ALL][i]->do_fsh_timestep_elem(time, u, elem, vCornerCoords);
244 : }
245 0 : UG_CATCH_THROW("DataEvaluatorBase::fsh_timestep_elem: Cannot finish timestep.");
246 0 : }
247 :
248 : template <typename TDomain>
249 0 : void DataEvaluator<TDomain>::
250 : add_jac_A_elem(LocalMatrix& J, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
251 : {
252 : UG_ASSERT(m_discPart & STIFF, "Using add_jac_A_elem, but not STIFF requested.");
253 :
254 : // compute elem-owned contribution
255 : try{
256 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
257 0 : m_vElemDisc[type][i]->do_add_jac_A_elem(J, u, elem, vCornerCoords);
258 : }
259 0 : UG_CATCH_THROW("DataEvaluatorBase::add_jac_A_elem: Cannot assemble Jacobian (A)");
260 :
261 : // compute linearized defect
262 : try{
263 0 : for(size_t i = 0; i < m_vImport[type][STIFF].size(); ++i)
264 0 : m_vImport[type][STIFF][i]->compute_lin_defect(u);
265 :
266 0 : for(size_t i = 0; i < m_vImport[type][RHS].size(); ++i)
267 0 : m_vImport[type][RHS][i]->compute_lin_defect(u);
268 : }
269 0 : UG_CATCH_THROW("DataEvaluatorBase::add_jac_A_elem: Cannot compute"
270 : " linearized defect for Import.");
271 :
272 : // add off diagonal coupling
273 : try{
274 : // loop all imports located in the stiffness part
275 0 : for(size_t i = 0; i < m_vImport[type][STIFF].size(); ++i)
276 0 : m_vImport[type][STIFF][i]->add_jacobian(J, 1.0);
277 :
278 : // loop all imports located in the rhs part
279 0 : for(size_t i = 0; i < m_vImport[type][RHS].size(); ++i)
280 0 : m_vImport[type][RHS][i]->add_jacobian(J, -1.0);
281 : }
282 0 : UG_CATCH_THROW("DataEvaluatorBase::add_jac_A_elem: Cannot add couplings.");
283 0 : }
284 :
285 : template <typename TDomain>
286 0 : void DataEvaluator<TDomain>::
287 : add_jac_M_elem(LocalMatrix& J, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
288 : {
289 : UG_ASSERT(m_discPart & MASS, "Using add_jac_M_elem, but not MASS requested.");
290 :
291 : // compute elem-owned contribution
292 : try{
293 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
294 0 : m_vElemDisc[type][i]->do_add_jac_M_elem(J, u, elem, vCornerCoords);
295 : }
296 0 : UG_CATCH_THROW("DataEvaluatorBase::add_jac_M_elem: Cannot assemble Jacobian (M)");
297 :
298 : // compute linearized defect
299 : try{
300 0 : for(size_t i = 0; i < m_vImport[type][MASS].size(); ++i)
301 0 : m_vImport[type][MASS][i]->compute_lin_defect(u);
302 : }
303 0 : UG_CATCH_THROW("DataEvaluatorBase::add_coupl_JM: Cannot compute"
304 : " linearized defect for Import.");
305 :
306 : // add off diagonal coupling
307 : try{
308 : // loop all imports located in the mass part
309 0 : for(size_t i = 0; i < m_vImport[type][MASS].size(); ++i)
310 0 : m_vImport[type][MASS][i]->add_jacobian(J, 1.0);
311 : }
312 0 : UG_CATCH_THROW("DataEvaluatorBase::add_coupl_JM: Cannot add couplings.");
313 0 : }
314 :
315 : template <typename TDomain>
316 0 : void DataEvaluator<TDomain>::
317 : add_def_A_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
318 : {
319 : UG_ASSERT(m_discPart & STIFF, "Using add_def_A_elem, but not STIFF requested.");
320 :
321 : try{
322 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
323 0 : m_vElemDisc[type][i]->do_add_def_A_elem(d, u, elem, vCornerCoords);
324 : }
325 0 : UG_CATCH_THROW("DataEvaluatorBase::add_def_A_elem: Cannot assemble Defect (A)");
326 0 : }
327 :
328 : template <typename TDomain>
329 0 : void DataEvaluator<TDomain>::
330 : add_def_A_expl_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
331 : {
332 : UG_ASSERT(m_discPart & EXPL, "Using add_def_A_elem, but not EXPL requested.");
333 :
334 : try{
335 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
336 0 : m_vElemDisc[type][i]->do_add_def_A_expl_elem(d, u, elem, vCornerCoords);
337 : }
338 0 : UG_CATCH_THROW("DataEvaluatorBase::add_def_A_expl_elem: Cannot assemble Defect (A)");
339 0 : }
340 :
341 : template <typename TDomain>
342 0 : void DataEvaluator<TDomain>::
343 : add_def_M_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
344 : {
345 : UG_ASSERT(m_discPart & MASS, "Using add_def_M_elem, but not MASS requested.");
346 :
347 : try{
348 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
349 0 : m_vElemDisc[type][i]->do_add_def_M_elem(d, u, elem, vCornerCoords);
350 : }
351 0 : UG_CATCH_THROW("DataEvaluatorBase::add_def_M_elem: Cannot assemble Defect (M)");
352 0 : }
353 :
354 : template <typename TDomain>
355 0 : void DataEvaluator<TDomain>::
356 : add_rhs_elem(LocalVector& rhs, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
357 : {
358 : UG_ASSERT(m_discPart & RHS, "Using add_rhs_elem, but not RHS requested.");
359 :
360 : try{
361 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
362 0 : m_vElemDisc[type][i]->do_add_rhs_elem(rhs, elem, vCornerCoords);
363 : }
364 0 : UG_CATCH_THROW("DataEvaluatorBase::add_rhs_elem: Cannot assemble rhs");
365 0 : }
366 :
367 : ////////////////////////////////////////////////////////////////////////////////
368 : // explicit template instantiations
369 : ////////////////////////////////////////////////////////////////////////////////
370 :
371 : #ifdef UG_DIM_1
372 : template class DataEvaluatorBase<Domain1d, IElemDisc<Domain1d> >;
373 : template class ErrorEvaluator<Domain1d>;
374 : template class DataEvaluator<Domain1d>;
375 : #endif
376 : #ifdef UG_DIM_2
377 : template class DataEvaluatorBase<Domain2d, IElemDisc<Domain2d> >;
378 : template class ErrorEvaluator<Domain2d>;
379 : template class DataEvaluator<Domain2d>;
380 : #endif
381 : #ifdef UG_DIM_3
382 : template class DataEvaluatorBase<Domain3d, IElemDisc<Domain3d> >;
383 : template class ErrorEvaluator<Domain3d>;
384 : template class DataEvaluator<Domain3d>;
385 : #endif
386 :
387 : } // end namespace ug
388 :
|
