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 "neumann_boundary_fv1.h"
34 : #include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
35 : #include "lib_disc/spatial_disc/disc_util/geom_provider.h"
36 :
37 : namespace ug{
38 :
39 :
40 : ////////////////////////////////////////////////////////////////////////////////
41 : // Constructor
42 : ////////////////////////////////////////////////////////////////////////////////
43 :
44 : template<typename TDomain>
45 0 : NeumannBoundaryFV1<TDomain>::NeumannBoundaryFV1(const char* function)
46 0 : :NeumannBoundaryBase<TDomain>(function)
47 : {
48 0 : register_all_funcs(false);
49 0 : }
50 :
51 : template<typename TDomain>
52 0 : void NeumannBoundaryFV1<TDomain>::
53 : prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid)
54 : {
55 0 : UG_COND_THROW(bNonRegularGrid && (TDomain::dim == 3),
56 : "NeumannBoundaryFV1: Hanging Nodes not implemented.");
57 :
58 0 : if(vLfeID.size() != 1)
59 0 : UG_THROW("NeumannBoundary: Need exactly 1 function.");
60 :
61 0 : if(vLfeID[0].order() != 1 || vLfeID[0].type() != LFEID::LAGRANGE)
62 0 : UG_THROW("NeumannBoundary: FV Scheme only implemented for 1st order Lagrange.");
63 0 : }
64 :
65 : template<typename TDomain>
66 0 : void NeumannBoundaryFV1<TDomain>::
67 : add(SmartPtr<CplUserData<number, dim> > data, const char* BndSubsets, const char* InnerSubsets)
68 : {
69 0 : m_vNumberData.push_back(NumberData(data, BndSubsets, InnerSubsets));
70 0 : this->add_inner_subsets(InnerSubsets);
71 0 : }
72 :
73 : template<typename TDomain>
74 0 : void NeumannBoundaryFV1<TDomain>::
75 : add(SmartPtr<CplUserData<number, dim, bool> > user, const char* BndSubsets, const char* InnerSubsets)
76 : {
77 0 : m_vBNDNumberData.push_back(BNDNumberData(user, BndSubsets, InnerSubsets));
78 0 : this->add_inner_subsets(InnerSubsets);
79 0 : }
80 :
81 : template<typename TDomain>
82 0 : void NeumannBoundaryFV1<TDomain>::
83 : add(SmartPtr<CplUserData<MathVector<dim>, dim> > user, const char* BndSubsets, const char* InnerSubsets)
84 : {
85 0 : m_vVectorData.push_back(VectorData(user, BndSubsets, InnerSubsets));
86 0 : this->add_inner_subsets(InnerSubsets);
87 0 : }
88 :
89 : template<typename TDomain>
90 0 : void NeumannBoundaryFV1<TDomain>::update_subset_groups()
91 : {
92 0 : for(size_t i = 0; i < m_vNumberData.size(); ++i)
93 0 : update_subset_groups(m_vNumberData[i]);
94 0 : for(size_t i = 0; i < m_vBNDNumberData.size(); ++i)
95 0 : update_subset_groups(m_vBNDNumberData[i]);
96 0 : for(size_t i = 0; i < m_vVectorData.size(); ++i)
97 0 : update_subset_groups(m_vVectorData[i]);
98 0 : }
99 :
100 : ////////////////////////////////////////////////////////////////////////////////
101 : // assembling functions
102 : ////////////////////////////////////////////////////////////////////////////////
103 :
104 : template<typename TDomain>
105 : template<typename TElem, typename TFVGeom>
106 0 : void NeumannBoundaryFV1<TDomain>::
107 : prep_elem_loop(const ReferenceObjectID roid, const int si)
108 : {
109 0 : update_subset_groups();
110 0 : m_si = si;
111 :
112 : // register subsetIndex at Geometry
113 0 : static TFVGeom& geo = GeomProvider<TFVGeom >::get();
114 :
115 : // request subset indices as boundary subset. This will force the
116 : // creation of boundary subsets when calling geo.update
117 :
118 0 : for(size_t i = 0; i < m_vNumberData.size(); ++i){
119 0 : if(!m_vNumberData[i].InnerSSGrp.contains(m_si)) continue;
120 0 : for(size_t s = 0; s < m_vNumberData[i].BndSSGrp.size(); ++s){
121 : const int si = m_vNumberData[i].BndSSGrp[s];
122 0 : geo.add_boundary_subset(si);
123 : }
124 : }
125 0 : for(size_t i = 0; i < m_vBNDNumberData.size(); ++i){
126 0 : if(!m_vBNDNumberData[i].InnerSSGrp.contains(m_si)) continue;
127 0 : for(size_t s = 0; s < m_vBNDNumberData[i].BndSSGrp.size(); ++s){
128 : const int si = m_vBNDNumberData[i].BndSSGrp[s];
129 0 : geo.add_boundary_subset(si);
130 : }
131 : }
132 0 : for(size_t i = 0; i < m_vVectorData.size(); ++i){
133 0 : if(!m_vVectorData[i].InnerSSGrp.contains(m_si)) continue;
134 0 : for(size_t s = 0; s < m_vVectorData[i].BndSSGrp.size(); ++s){
135 : const int si = m_vVectorData[i].BndSSGrp[s];
136 0 : geo.add_boundary_subset(si);
137 : }
138 : }
139 :
140 : // clear imports, since we will set them afterwards
141 : this->clear_imports();
142 :
143 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
144 :
145 : // set lin defect fct for imports
146 0 : for(size_t data = 0; data < m_vNumberData.size(); ++data)
147 : {
148 0 : if(!m_vNumberData[data].InnerSSGrp.contains(m_si)) continue;
149 0 : m_vNumberData[data].import.set_fct(id,
150 : &m_vNumberData[data],
151 : &NumberData::template lin_def<TElem, TFVGeom>);
152 :
153 0 : this->register_import(m_vNumberData[data].import);
154 : m_vNumberData[data].import.set_rhs_part();
155 : }
156 0 : }
157 :
158 : template<typename TDomain>
159 : template<typename TElem, typename TFVGeom>
160 0 : void NeumannBoundaryFV1<TDomain>::
161 : prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[])
162 : {
163 : // update Geometry for this element
164 0 : static TFVGeom& geo = GeomProvider<TFVGeom >::get();
165 : try{
166 0 : geo.update(elem, vCornerCoords, &(this->subset_handler()));
167 : }
168 0 : UG_CATCH_THROW("NeumannBoundaryFV1::prep_elem: "
169 : "Cannot update Finite Volume Geometry.");
170 :
171 0 : for(size_t i = 0; i < m_vNumberData.size(); ++i)
172 0 : if(m_vNumberData[i].InnerSSGrp.contains(m_si))
173 0 : m_vNumberData[i].template extract_bip<TElem, TFVGeom>(geo);
174 0 : }
175 :
176 : template<typename TDomain>
177 : template<typename TElem, typename TFVGeom>
178 0 : void NeumannBoundaryFV1<TDomain>::
179 : add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[])
180 : {
181 0 : const static TFVGeom& geo = GeomProvider<TFVGeom >::get();
182 : typedef typename TFVGeom::BF BF;
183 :
184 : // Number Data
185 0 : for(size_t data = 0; data < m_vNumberData.size(); ++data){
186 0 : if(!m_vNumberData[data].InnerSSGrp.contains(m_si)) continue;
187 : size_t ip = 0;
188 0 : for(size_t s = 0; s < m_vNumberData[data].BndSSGrp.size(); ++s){
189 : const int si = m_vNumberData[data].BndSSGrp[s];
190 0 : const std::vector<BF>& vBF = geo.bf(si);
191 :
192 0 : for(size_t i = 0; i < vBF.size(); ++i, ++ip){
193 0 : const int co = vBF[i].node_id();
194 0 : d(_C_, co) -= m_vNumberData[data].import[ip]
195 0 : * vBF[i].volume();
196 : }
197 : }
198 : }
199 :
200 : // conditional Number Data
201 0 : for(size_t data = 0; data < m_vBNDNumberData.size(); ++data){
202 0 : if(!m_vBNDNumberData[data].InnerSSGrp.contains(m_si)) continue;
203 0 : for(size_t s = 0; s < m_vBNDNumberData[data].BndSSGrp.size(); ++s) {
204 : const int si = m_vBNDNumberData[data].BndSSGrp[s];
205 0 : const std::vector<BF>& vBF = geo.bf(si);
206 :
207 0 : for(size_t i = 0; i < vBF.size(); ++i){
208 0 : number val = 0.0;
209 0 : if(!(*m_vBNDNumberData[data].functor)(val, vBF[i].global_ip(), this->time(), si))
210 0 : continue;
211 :
212 0 : const int co = vBF[i].node_id();
213 0 : d(_C_, co) -= val * vBF[i].volume();
214 : }
215 : }
216 : }
217 :
218 : // vector data
219 0 : for(size_t data = 0; data < m_vVectorData.size(); ++data){
220 0 : if(!m_vVectorData[data].InnerSSGrp.contains(m_si)) continue;
221 0 : for(size_t s = 0; s < m_vVectorData[data].BndSSGrp.size(); ++s){
222 : const int si = m_vVectorData[data].BndSSGrp[s];
223 0 : const std::vector<BF>& vBF = geo.bf(si);
224 :
225 0 : for(size_t i = 0; i < vBF.size(); ++i){
226 : MathVector<dim> val;
227 0 : (*m_vVectorData[data].functor)(val, vBF[i].global_ip(), this->time(), si);
228 :
229 0 : const int co = vBF[i].node_id();
230 0 : d(_C_, co) -= VecDot(val, vBF[i].normal());
231 : }
232 : }
233 : }
234 0 : }
235 :
236 : template<typename TDomain>
237 : template<typename TElem, typename TGeom>
238 0 : void NeumannBoundaryFV1<TDomain>::
239 : fsh_elem_loop()
240 : {
241 : // remove subsetIndex from Geometry
242 0 : static TGeom& geo = GeomProvider<TGeom >::get();
243 :
244 :
245 : // unrequest subset indices as boundary subset. This will force the
246 : // creation of boundary subsets when calling geo.update
247 :
248 0 : for(size_t i = 0; i < m_vNumberData.size(); ++i){
249 0 : if(!m_vNumberData[i].InnerSSGrp.contains(m_si)) continue;
250 0 : for(size_t s = 0; s < m_vNumberData[i].BndSSGrp.size(); ++s){
251 : const int si = m_vNumberData[i].BndSSGrp[s];
252 0 : geo.remove_boundary_subset(si);
253 0 : geo.reset_curr_elem();
254 : }
255 : }
256 :
257 0 : for(size_t i = 0; i < m_vBNDNumberData.size(); ++i){
258 0 : if(!m_vBNDNumberData[i].InnerSSGrp.contains(m_si)) continue;
259 0 : for(size_t s = 0; s < m_vBNDNumberData[i].BndSSGrp.size(); ++s){
260 : const int si = m_vBNDNumberData[i].BndSSGrp[s];
261 0 : geo.remove_boundary_subset(si);
262 0 : geo.reset_curr_elem();
263 : }
264 : }
265 :
266 0 : for(size_t i = 0; i < m_vVectorData.size(); ++i){
267 0 : if(!m_vVectorData[i].InnerSSGrp.contains(m_si)) continue;
268 0 : for(size_t s = 0; s < m_vVectorData[i].BndSSGrp.size(); ++s){
269 : const int si = m_vVectorData[i].BndSSGrp[s];
270 0 : geo.remove_boundary_subset(si);
271 0 : geo.reset_curr_elem();
272 : }
273 : }
274 0 : }
275 :
276 :
277 : ///////////////////////////////////////////////////////////////////////////////////////////////////
278 : // error estimation (begin) //
279 :
280 : // prepares the loop over all elements of one type for the computation of the error estimator
281 : template<typename TDomain>
282 : template <typename TElem, typename TFVGeom>
283 0 : void NeumannBoundaryFV1<TDomain>::
284 : prep_err_est_elem_loop(const ReferenceObjectID roid, const int si)
285 : {
286 : // get the error estimator data object and check that it is of the right type
287 : // we check this at this point in order to be able to dispense with this check later on
288 : // (i.e. in prep_err_est_elem and compute_err_est_A_elem())
289 0 : if (this->m_spErrEstData.get() == NULL)
290 : {
291 0 : UG_THROW("No ErrEstData object has been given to this ElemDisc!");
292 : }
293 :
294 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
295 :
296 0 : if (!err_est_data)
297 : {
298 0 : UG_THROW("Dynamic cast to SideAndElemErrEstData failed."
299 : << std::endl << "Make sure you handed the correct type of ErrEstData to this discretization.");
300 : }
301 :
302 0 : update_subset_groups();
303 0 : m_si = si;
304 :
305 : // clear imports, since we will set them now
306 : this->clear_imports();
307 :
308 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
309 :
310 : // set lin defect fct for imports
311 0 : for (size_t data = 0; data < m_vNumberData.size(); ++data)
312 : {
313 0 : if (!m_vNumberData[data].InnerSSGrp.contains(m_si)) continue;
314 0 : m_vNumberData[data].import.set_fct(id,
315 : &m_vNumberData[data],
316 : &NumberData::template lin_def<TElem, TFVGeom>);
317 :
318 0 : this->register_import(m_vNumberData[data].import);
319 : m_vNumberData[data].import.set_rhs_part();
320 : }
321 :
322 : // get local IPs for imports
323 : static const int refDim = TElem::dim;
324 :
325 : size_t numSideIPs;
326 : const MathVector<refDim>* sideIPs;
327 : try
328 : {
329 : numSideIPs = err_est_data->num_all_side_ips(roid);
330 0 : sideIPs = err_est_data->template side_local_ips<refDim>(roid);
331 :
332 0 : if (!sideIPs) return; // is NULL if TElem is not of the same dim as TDomain
333 : }
334 0 : UG_CATCH_THROW("Integration points for error estimator cannot be set.");
335 :
336 0 : for (size_t i = 0; i < m_vNumberData.size(); ++i)
337 : {
338 0 : if (m_vNumberData[i].InnerSSGrp.contains(m_si))
339 : m_vNumberData[i].template set_local_ips<refDim>(sideIPs, numSideIPs);
340 : }
341 : };
342 :
343 : template<typename TDomain>
344 : template<typename TElem, typename TFVGeom>
345 0 : void NeumannBoundaryFV1<TDomain>::
346 : prep_err_est_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
347 : {
348 : // get global IPs for imports
349 : size_t numSideIPs;
350 : MathVector<dim>* sideIPs;
351 :
352 : try
353 : {
354 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
355 :
356 0 : numSideIPs = err_est_data->num_all_side_ips(elem->reference_object_id());
357 0 : sideIPs = err_est_data->all_side_global_ips(elem, vCornerCoords);
358 : }
359 0 : UG_CATCH_THROW("Global integration points for error estimator cannot be set.");
360 :
361 :
362 0 : for (size_t i = 0; i < m_vNumberData.size(); ++i)
363 : {
364 0 : if (m_vNumberData[i].InnerSSGrp.contains(m_si))
365 : m_vNumberData[i].set_global_ips(&sideIPs[0], numSideIPs);
366 : }
367 0 : }
368 :
369 : // computes the error estimator contribution for one element
370 : template<typename TDomain>
371 : template <typename TElem, typename TFVGeom>
372 0 : void NeumannBoundaryFV1<TDomain>::
373 : compute_err_est_rhs_elem(GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale)
374 : {
375 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
376 :
377 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
378 :
379 0 : if (err_est_data->surface_view().get() == NULL) {UG_THROW("Error estimator has NULL surface view.");}
380 0 : MultiGrid* pErrEstGrid = (MultiGrid*) (err_est_data->surface_view()->subset_handler()->multi_grid());
381 :
382 : // get the sides of the element
383 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::side_type>::secure_container side_list;
384 0 : pErrEstGrid->associated_elements_sorted(side_list, (TElem*) elem);
385 0 : if (side_list.size() != (size_t) ref_elem_type::numSides)
386 0 : UG_THROW ("Mismatch of numbers of sides in 'ConvectionDiffusionFV1::compute_err_est_elem'");
387 :
388 : // Number Data
389 0 : for (size_t data = 0; data < m_vNumberData.size(); ++data)
390 : {
391 0 : if (!m_vNumberData[data].InnerSSGrp.contains(m_si)) continue;
392 :
393 : // loop sides
394 0 : for (size_t side = 0; side < side_list.size(); side++)
395 : {
396 0 : for (size_t ss = 0; ss < m_vNumberData[data].BndSSGrp.size(); ss++)
397 : {
398 : // is the bnd cond subset index the same as the current side's?
399 0 : if (err_est_data->surface_view()->subset_handler()->get_subset_index(side_list[side])
400 0 : != m_vNumberData[data].BndSSGrp[ss]) continue;
401 :
402 0 : const size_t ipIndexStart = err_est_data->first_side_ips(elem->reference_object_id(), side);
403 :
404 0 : for (size_t sip = 0; sip < err_est_data->num_side_ips(side_list[side]); sip++)
405 : {
406 0 : size_t ip = ipIndexStart + sip;
407 : // sign must be negative, since the data represents efflux
408 0 : (*err_est_data)(side_list[side],sip) -= scale * m_vNumberData[data].import[ip];
409 : }
410 : }
411 : }
412 : }
413 :
414 : // store global side IPs
415 0 : ReferenceObjectID roid = elem->reference_object_id();
416 0 : MathVector<dim>* glob_ips = err_est_data->all_side_global_ips(elem, vCornerCoords);
417 :
418 : // conditional Number Data
419 0 : for (size_t data = 0; data < m_vBNDNumberData.size(); ++data)
420 : {
421 0 : if (!m_vBNDNumberData[data].InnerSSGrp.contains(m_si)) continue;
422 :
423 : // loop sides
424 0 : for (size_t side = 0; side < side_list.size(); side++)
425 : {
426 0 : for (size_t ss = 0; ss < m_vBNDNumberData[data].BndSSGrp.size(); ss++)
427 : {
428 : // is the bnd cond subset index the same as the current side's?
429 0 : if (err_est_data->surface_view()->subset_handler()->get_subset_index(side_list[side])
430 0 : != m_vBNDNumberData[data].BndSSGrp[ss]) continue;
431 :
432 : int si = m_vBNDNumberData[data].BndSSGrp[ss];
433 :
434 : const size_t ipIndexStart = err_est_data->first_side_ips(roid, side);
435 :
436 0 : for (size_t sip = 0; sip < err_est_data->num_side_ips(side_list[side]); sip++)
437 : {
438 0 : size_t ip = ipIndexStart + sip;
439 0 : number val = 0.0;
440 0 : if (!(*m_vBNDNumberData[data].functor)(val, glob_ips[ip], this->time(), si)) continue;
441 :
442 : // sign must be negative, since the data represents efflux
443 0 : (*err_est_data)(side_list[side],sip) -= scale * val;
444 : }
445 : }
446 : }
447 : }
448 :
449 : // vector data
450 0 : for (size_t data = 0; data < m_vVectorData.size(); ++data)
451 : {
452 0 : if (!m_vVectorData[data].InnerSSGrp.contains(m_si)) continue;
453 :
454 : // loop sides
455 0 : for (size_t side = 0; side < side_list.size(); side++)
456 : {
457 : MathVector<dim> fd, normal;
458 :
459 : // normal on side
460 0 : SideNormal<ref_elem_type,dim>(normal, side, vCornerCoords);
461 0 : VecNormalize(normal, normal);
462 :
463 0 : for (size_t ss = 0; ss < m_vVectorData[data].BndSSGrp.size(); ss++)
464 : {
465 : // is the bnd cond subset index the same as the current side's?
466 0 : if (err_est_data->surface_view()->subset_handler()->get_subset_index(side_list[side])
467 0 : != m_vVectorData[data].BndSSGrp[ss]) continue;
468 :
469 : int si = m_vBNDNumberData[data].BndSSGrp[ss];
470 :
471 : const size_t ipIndexStart = err_est_data->first_side_ips(roid, side);
472 :
473 0 : for (size_t sip = 0; sip < err_est_data->num_side_ips(side_list[side]); sip++)
474 : {
475 0 : size_t ip = ipIndexStart + sip;
476 0 : (*m_vVectorData[data].functor)(fd, glob_ips[ip], this->time(), si);
477 :
478 : // sign must be negative, since the data represents efflux
479 0 : (*err_est_data)(side_list[side],sip) -= scale * VecDot(fd, normal);
480 : }
481 : }
482 : }
483 : }
484 0 : }
485 :
486 : // postprocesses the loop over all elements of one type in the computation of the error estimator
487 : template<typename TDomain>
488 : template <typename TElem, typename TFVGeom>
489 0 : void NeumannBoundaryFV1<TDomain>::
490 : fsh_err_est_elem_loop()
491 : {
492 : // finish the element loop in the same way as the actual discretization
493 0 : this->template fsh_elem_loop<TElem, TFVGeom> ();
494 0 : }
495 :
496 : // error estimation (end) //
497 : ///////////////////////////////////////////////////////////////////////////////////////////////////
498 :
499 :
500 : ////////////////////////////////////////////////////////////////////////////////
501 : // Number Data
502 : ////////////////////////////////////////////////////////////////////////////////
503 :
504 :
505 : template<typename TDomain>
506 : template<typename TElem, typename TFVGeom>
507 0 : void NeumannBoundaryFV1<TDomain>::NumberData::
508 : lin_def(const LocalVector& u,
509 : std::vector<std::vector<number> > vvvLinDef[],
510 : const size_t nip)
511 : {
512 : // get finite volume geometry
513 0 : const static TFVGeom& geo = GeomProvider<TFVGeom>::get();
514 : typedef typename TFVGeom::BF BF;
515 :
516 0 : for(size_t s = 0; s < this->BndSSGrp.size(); ++s)
517 : {
518 : const int si = this->BndSSGrp[s];
519 0 : const std::vector<BF>& vBF = geo.bf(si);
520 0 : for(size_t i = 0; i < vBF.size(); ++i){
521 0 : const int co = vBF[i].node_id();
522 0 : vvvLinDef[i][_C_][co] -= vBF[i].volume();
523 : }
524 : }
525 0 : }
526 :
527 : template<typename TDomain>
528 : template<typename TElem, typename TFVGeom>
529 0 : void NeumannBoundaryFV1<TDomain>::NumberData::
530 : extract_bip(const TFVGeom& geo)
531 : {
532 : typedef typename TFVGeom::BF BF;
533 : static const int locDim = TElem::dim;
534 :
535 : std::vector<MathVector<locDim> >* vLocIP = local_ips<locDim>();
536 :
537 : vLocIP->clear();
538 : vGloIP.clear();
539 0 : for(size_t s = 0; s < this->BndSSGrp.size(); s++)
540 : {
541 : const int si = this->BndSSGrp[s];
542 : const std::vector<BF>& vBF = geo.bf(si);
543 0 : for(size_t i = 0; i < vBF.size(); ++i)
544 : {
545 : const BF& bf = vBF[i];
546 0 : vLocIP->push_back(bf.local_ip());
547 0 : vGloIP.push_back(bf.global_ip());
548 : }
549 : }
550 :
551 : // either both are empty or none is empty
552 : UG_ASSERT((vLocIP->empty() && vGloIP.empty()) || (!(vLocIP->empty() || vGloIP.empty())),
553 : "Either vLocIP and vGloIP both have to be empty or both have to be filled!");
554 :
555 0 : if(vLocIP->empty()){
556 0 : import.template set_local_ips<locDim>(NULL, 0);
557 0 : import.set_global_ips(NULL, 0);
558 : }
559 : else{
560 0 : import.template set_local_ips<locDim>(&(*vLocIP)[0], vLocIP->size());
561 0 : import.set_global_ips(&vGloIP[0], vGloIP.size());
562 : }
563 0 : }
564 :
565 : template<typename TDomain>
566 : template <int refDim>
567 : std::vector<MathVector<refDim> >* NeumannBoundaryFV1<TDomain>::NumberData::local_ips()
568 : {
569 : switch (refDim)
570 : {
571 0 : case 1: return (std::vector<MathVector<refDim> >*)(&vLocIP_dim1);
572 0 : case 2: return (std::vector<MathVector<refDim> >*)(&vLocIP_dim2);
573 0 : case 3: return (std::vector<MathVector<refDim> >*)(&vLocIP_dim3);
574 : }
575 : }
576 :
577 : ////////////////////////////////////////////////////////////////////////////////
578 : // register assemble functions
579 : ////////////////////////////////////////////////////////////////////////////////
580 :
581 : #ifdef UG_DIM_1
582 : template<>
583 0 : void NeumannBoundaryFV1<Domain1d>::register_all_funcs(bool bHang)
584 : {
585 0 : register_func<RegularEdge, FV1Geometry<RegularEdge, dim> >();
586 0 : }
587 : #endif
588 :
589 : #ifdef UG_DIM_2
590 : template<>
591 0 : void NeumannBoundaryFV1<Domain2d>::register_all_funcs(bool bHang)
592 : {
593 0 : register_func<RegularEdge, FV1Geometry<RegularEdge, dim> >();
594 0 : register_func<Triangle, FV1Geometry<Triangle, dim> >();
595 0 : register_func<Quadrilateral, FV1Geometry<Quadrilateral, dim> >();
596 0 : }
597 : #endif
598 :
599 : #ifdef UG_DIM_3
600 : template<>
601 0 : void NeumannBoundaryFV1<Domain3d>::register_all_funcs(bool bHang)
602 : {
603 0 : if(!bHang){
604 0 : register_func<RegularEdge, FV1Geometry<RegularEdge, dim> >();
605 0 : register_func<Triangle, FV1Geometry<Triangle, dim> >();
606 0 : register_func<Quadrilateral, FV1Geometry<Quadrilateral, dim> >();
607 0 : register_func<Tetrahedron, FV1Geometry<Tetrahedron, dim> >();
608 0 : register_func<Prism, FV1Geometry<Prism, dim> >();
609 0 : register_func<Pyramid, FV1Geometry<Pyramid, dim> >();
610 0 : register_func<Hexahedron, FV1Geometry<Hexahedron, dim> >();
611 0 : register_func<Octahedron, FV1Geometry<Octahedron, dim> >();
612 : }
613 0 : else UG_THROW("NeumannBoundaryFV1: Hanging Nodes not implemented.");
614 0 : }
615 : #endif
616 :
617 : template<typename TDomain>
618 : template<typename TElem, typename TFVGeom>
619 0 : void NeumannBoundaryFV1<TDomain>::register_func()
620 : {
621 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
622 : typedef this_type T;
623 :
624 : this->clear_add_fct(id);
625 :
626 : this->set_prep_elem_loop_fct(id, &T::template prep_elem_loop<TElem, TFVGeom>);
627 : this->set_prep_elem_fct( id, &T::template prep_elem<TElem, TFVGeom>);
628 : this->set_add_rhs_elem_fct( id, &T::template add_rhs_elem<TElem, TFVGeom>);
629 : this->set_fsh_elem_loop_fct( id, &T::template fsh_elem_loop<TElem, TFVGeom>);
630 :
631 : this->set_add_jac_A_elem_fct( id, &T::template add_jac_A_elem<TElem, TFVGeom>);
632 : this->set_add_jac_M_elem_fct( id, &T::template add_jac_M_elem<TElem, TFVGeom>);
633 : this->set_add_def_A_elem_fct( id, &T::template add_def_A_elem<TElem, TFVGeom>);
634 : this->set_add_def_M_elem_fct( id, &T::template add_def_M_elem<TElem, TFVGeom>);
635 :
636 : // error estimator parts
637 : this->set_prep_err_est_elem_loop(id, &T::template prep_err_est_elem_loop<TElem, TFVGeom>);
638 : this->set_prep_err_est_elem(id, &T::template prep_err_est_elem<TElem, TFVGeom>);
639 : this->set_compute_err_est_rhs_elem(id, &T::template compute_err_est_rhs_elem<TElem, TFVGeom>);
640 : this->set_fsh_err_est_elem_loop(id, &T::template fsh_err_est_elem_loop<TElem, TFVGeom>);
641 0 : }
642 :
643 : ////////////////////////////////////////////////////////////////////////////////
644 : // explicit template instantiations
645 : ////////////////////////////////////////////////////////////////////////////////
646 :
647 : #ifdef UG_DIM_1
648 : template class NeumannBoundaryFV1<Domain1d>;
649 : #endif
650 : #ifdef UG_DIM_2
651 : template class NeumannBoundaryFV1<Domain2d>;
652 : #endif
653 : #ifdef UG_DIM_3
654 : template class NeumannBoundaryFV1<Domain3d>;
655 : #endif
656 :
657 : } // namespace ug
658 :
|
