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 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i){
166 : 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);
167 0 : m_vElemDisc[PT_ALL][i]->do_prep_elem(u, elem, roid, vCornerCoords);
168 : }
169 : }
170 0 : UG_CATCH_THROW("DataEvaluatorBase::prep_elem: Cannot prepare element.");
171 :
172 : // adjust lin defect array of imports and derivative array of exports
173 : // INFO: This is place here, since the 'prepare_elem' method of an element
174 : // disc may change the number of integration points, even if the type
175 : // of the element (e.g. triangle, quad) stays the same. This is the
176 : // case for, e.g., the NeumannBoundary element disc.
177 0 : if(bDeriv)
178 : {
179 0 : for(size_t i = 0; i < m_vImport[PT_ALL][MASS].size(); ++i)
180 0 : m_vImport[PT_ALL][MASS][i]->update_dof_sizes(ind);
181 0 : for(size_t i = 0; i < m_vImport[PT_ALL][STIFF].size(); ++i)
182 0 : m_vImport[PT_ALL][STIFF][i]->update_dof_sizes(ind);
183 0 : for(size_t i = 0; i < m_vImport[PT_ALL][RHS].size(); ++i)
184 0 : m_vImport[PT_ALL][RHS][i]->update_dof_sizes(ind);
185 :
186 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i)
187 0 : m_vDependentData[i]->update_dof_sizes(ind);
188 : }
189 :
190 : // evaluate position data
191 0 : for(size_t i = 0; i < m_vPosData.size(); ++i)
192 0 : m_vPosData[i]->compute(&u, elem, vCornerCoords, false);
193 :
194 : // process dependent data:
195 : // We can not simply compute exports first, then Linker, because an export
196 : // itself could depend on other data if implemented somehow in the IElemDisc
197 : // (e.g. using data from some DataImport). Thus, we have to loop the sorted
198 : // vector of all dependent data (that is correctly sorted the way that always
199 : // needed data has previously computed).
200 :
201 : // compute the data
202 : try{
203 0 : if (! time_series_needed ()) { // assemble for the given LocalVector
204 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i){
205 0 : u.access_by_map(m_vDependentData[i]->map());
206 0 : m_vDependentData[i]->compute(&u, elem, vCornerCoords, bDeriv);
207 : }
208 : }
209 : else { // assemble for LocalVectorTimeSeries
210 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i){
211 0 : u.access_by_map(m_vDependentData[i]->map());
212 0 : m_vDependentData[i]->compute(m_pLocTimeSeries, elem, vCornerCoords, bDeriv);
213 : }
214 : }
215 : }
216 0 : UG_CATCH_THROW("DataEvaluatorBase::prep_elem: Cannot compute data for Export or Linker.");
217 0 : }
218 :
219 : template <typename TDomain>
220 0 : void DataEvaluator<TDomain>::finish_timestep(const number time, VectorProxyBase* u, size_t algebra_id)
221 : {
222 : try
223 : {
224 0 : for (size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
225 0 : m_vElemDisc[PT_ALL][i]->do_fsh_timestep(time, u, algebra_id);
226 : }
227 0 : UG_CATCH_THROW("DataEvaluatorBase::finish_timestep: Cannot prepare time step.");
228 0 : }
229 :
230 : template <typename TDomain>
231 0 : void DataEvaluator<TDomain>::
232 : finish_timestep_elem(const number time, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
233 : {
234 : try{
235 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
236 0 : m_vElemDisc[PT_ALL][i]->do_fsh_timestep_elem(time, u, elem, vCornerCoords);
237 : }
238 0 : UG_CATCH_THROW("DataEvaluatorBase::fsh_timestep_elem: Cannot finish timestep.");
239 0 : }
240 :
241 : template <typename TDomain>
242 0 : void DataEvaluator<TDomain>::
243 : add_jac_A_elem(LocalMatrix& J, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
244 : {
245 : UG_ASSERT(m_discPart & STIFF, "Using add_jac_A_elem, but not STIFF requested.");
246 :
247 : // compute elem-owned contribution
248 : try{
249 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
250 0 : m_vElemDisc[type][i]->do_add_jac_A_elem(J, u, elem, vCornerCoords);
251 : }
252 0 : UG_CATCH_THROW("DataEvaluatorBase::add_jac_A_elem: Cannot assemble Jacobian (A)");
253 :
254 : // compute linearized defect
255 : try{
256 0 : for(size_t i = 0; i < m_vImport[type][STIFF].size(); ++i)
257 0 : m_vImport[type][STIFF][i]->compute_lin_defect(u);
258 :
259 0 : for(size_t i = 0; i < m_vImport[type][RHS].size(); ++i)
260 0 : m_vImport[type][RHS][i]->compute_lin_defect(u);
261 : }
262 0 : UG_CATCH_THROW("DataEvaluatorBase::add_jac_A_elem: Cannot compute"
263 : " linearized defect for Import.");
264 :
265 : // add off diagonal coupling
266 : try{
267 : // loop all imports located in the stiffness part
268 0 : for(size_t i = 0; i < m_vImport[type][STIFF].size(); ++i)
269 0 : m_vImport[type][STIFF][i]->add_jacobian(J, 1.0);
270 :
271 : // loop all imports located in the rhs part
272 0 : for(size_t i = 0; i < m_vImport[type][RHS].size(); ++i)
273 0 : m_vImport[type][RHS][i]->add_jacobian(J, -1.0);
274 : }
275 0 : UG_CATCH_THROW("DataEvaluatorBase::add_jac_A_elem: Cannot add couplings.");
276 0 : }
277 :
278 : template <typename TDomain>
279 0 : void DataEvaluator<TDomain>::
280 : add_jac_M_elem(LocalMatrix& J, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
281 : {
282 : UG_ASSERT(m_discPart & MASS, "Using add_jac_M_elem, but not MASS requested.");
283 :
284 : // compute elem-owned contribution
285 : try{
286 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
287 0 : m_vElemDisc[type][i]->do_add_jac_M_elem(J, u, elem, vCornerCoords);
288 : }
289 0 : UG_CATCH_THROW("DataEvaluatorBase::add_jac_M_elem: Cannot assemble Jacobian (M)");
290 :
291 : // compute linearized defect
292 : try{
293 0 : for(size_t i = 0; i < m_vImport[type][MASS].size(); ++i)
294 0 : m_vImport[type][MASS][i]->compute_lin_defect(u);
295 : }
296 0 : UG_CATCH_THROW("DataEvaluatorBase::add_coupl_JM: Cannot compute"
297 : " linearized defect for Import.");
298 :
299 : // add off diagonal coupling
300 : try{
301 : // loop all imports located in the mass part
302 0 : for(size_t i = 0; i < m_vImport[type][MASS].size(); ++i)
303 0 : m_vImport[type][MASS][i]->add_jacobian(J, 1.0);
304 : }
305 0 : UG_CATCH_THROW("DataEvaluatorBase::add_coupl_JM: Cannot add couplings.");
306 0 : }
307 :
308 : template <typename TDomain>
309 0 : void DataEvaluator<TDomain>::
310 : add_def_A_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
311 : {
312 : UG_ASSERT(m_discPart & STIFF, "Using add_def_A_elem, but not STIFF requested.");
313 :
314 : try{
315 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
316 0 : m_vElemDisc[type][i]->do_add_def_A_elem(d, u, elem, vCornerCoords);
317 : }
318 0 : UG_CATCH_THROW("DataEvaluatorBase::add_def_A_elem: Cannot assemble Defect (A)");
319 0 : }
320 :
321 : template <typename TDomain>
322 0 : void DataEvaluator<TDomain>::
323 : add_def_A_expl_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
324 : {
325 : UG_ASSERT(m_discPart & EXPL, "Using add_def_A_elem, but not EXPL requested.");
326 :
327 : try{
328 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
329 0 : m_vElemDisc[type][i]->do_add_def_A_expl_elem(d, u, elem, vCornerCoords);
330 : }
331 0 : UG_CATCH_THROW("DataEvaluatorBase::add_def_A_expl_elem: Cannot assemble Defect (A)");
332 0 : }
333 :
334 : template <typename TDomain>
335 0 : void DataEvaluator<TDomain>::
336 : add_def_M_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
337 : {
338 : UG_ASSERT(m_discPart & MASS, "Using add_def_M_elem, but not MASS requested.");
339 :
340 : try{
341 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
342 0 : m_vElemDisc[type][i]->do_add_def_M_elem(d, u, elem, vCornerCoords);
343 : }
344 0 : UG_CATCH_THROW("DataEvaluatorBase::add_def_M_elem: Cannot assemble Defect (M)");
345 0 : }
346 :
347 : template <typename TDomain>
348 0 : void DataEvaluator<TDomain>::
349 : add_rhs_elem(LocalVector& rhs, GridObject* elem, const MathVector<dim> vCornerCoords[], ProcessType type)
350 : {
351 : UG_ASSERT(m_discPart & RHS, "Using add_rhs_elem, but not RHS requested.");
352 :
353 : try{
354 0 : for(size_t i = 0; i < m_vElemDisc[type].size(); ++i)
355 0 : m_vElemDisc[type][i]->do_add_rhs_elem(rhs, elem, vCornerCoords);
356 : }
357 0 : UG_CATCH_THROW("DataEvaluatorBase::add_rhs_elem: Cannot assemble rhs");
358 0 : }
359 :
360 : ////////////////////////////////////////////////////////////////////////////////
361 : // explicit template instantiations
362 : ////////////////////////////////////////////////////////////////////////////////
363 :
364 : #ifdef UG_DIM_1
365 : template class DataEvaluatorBase<Domain1d, IElemDisc<Domain1d> >;
366 : template class ErrorEvaluator<Domain1d>;
367 : template class DataEvaluator<Domain1d>;
368 : #endif
369 : #ifdef UG_DIM_2
370 : template class DataEvaluatorBase<Domain2d, IElemDisc<Domain2d> >;
371 : template class ErrorEvaluator<Domain2d>;
372 : template class DataEvaluator<Domain2d>;
373 : #endif
374 : #ifdef UG_DIM_3
375 : template class DataEvaluatorBase<Domain3d, IElemDisc<Domain3d> >;
376 : template class ErrorEvaluator<Domain3d>;
377 : template class DataEvaluator<Domain3d>;
378 : #endif
379 :
380 : } // end namespace ug
381 :
|
