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_fv.h"
34 : #include "lib_disc/spatial_disc/disc_util/fvho_geom.h"
35 : #include "lib_disc/spatial_disc/disc_util/geom_provider.h"
36 :
37 : namespace ug{
38 :
39 : ////////////////////////////////////////////////////////////////////////////////
40 : // Constructor
41 : ////////////////////////////////////////////////////////////////////////////////
42 :
43 : template<typename TDomain>
44 0 : NeumannBoundaryFV<TDomain>::NeumannBoundaryFV(const char* function)
45 : :NeumannBoundaryBase<TDomain>(function),
46 0 : m_order(1), m_lfeID(LFEID::LAGRANGE, TDomain::dim, m_order)
47 : {
48 0 : register_all_funcs(m_order);
49 0 : }
50 :
51 : template<typename TDomain>
52 0 : void NeumannBoundaryFV<TDomain>::
53 : prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid)
54 : {
55 0 : if(bNonRegularGrid)
56 0 : UG_THROW("NeumannBoundary: Hanging Nodes not implemented.");
57 :
58 : // check number
59 0 : if(vLfeID.size() != 1)
60 0 : UG_THROW("NeumannBoundaryFV: Need exactly 1 function.");
61 :
62 0 : if(vLfeID[0].type() != LFEID::LAGRANGE)
63 0 : UG_THROW("NeumannBoundary: FV Scheme only implemented for Lagrange-type elements.");
64 :
65 : // check that not ADAPTIVE
66 0 : if(vLfeID[0].order() < 1)
67 0 : UG_THROW("NeumannBoundary: Adaptive order not implemented.");
68 :
69 : // set order
70 0 : m_lfeID = vLfeID[0];
71 0 : m_order = vLfeID[0].order();
72 :
73 0 : register_all_funcs(m_order);
74 0 : }
75 :
76 : template<typename TDomain>
77 0 : void NeumannBoundaryFV<TDomain>::
78 : add(SmartPtr<CplUserData<number, dim> > data, const char* BndSubsets, const char* InnerSubsets)
79 : {
80 0 : m_vNumberData.push_back(NumberData(data, BndSubsets, InnerSubsets, this));
81 0 : this->add_inner_subsets(InnerSubsets);
82 0 : }
83 :
84 : template<typename TDomain>
85 0 : void NeumannBoundaryFV<TDomain>::
86 : add(SmartPtr<CplUserData<number, dim, bool> > user, const char* BndSubsets, const char* InnerSubsets)
87 : {
88 0 : m_vBNDNumberData.push_back(BNDNumberData(user, BndSubsets, InnerSubsets));
89 0 : this->add_inner_subsets(InnerSubsets);
90 0 : }
91 :
92 : template<typename TDomain>
93 0 : void NeumannBoundaryFV<TDomain>::
94 : add(SmartPtr<CplUserData<MathVector<dim>, dim> > user, const char* BndSubsets, const char* InnerSubsets)
95 : {
96 0 : m_vVectorData.push_back(VectorData(user, BndSubsets, InnerSubsets));
97 0 : this->add_inner_subsets(InnerSubsets);
98 0 : }
99 :
100 : template<typename TDomain>
101 0 : void NeumannBoundaryFV<TDomain>::update_subset_groups()
102 : {
103 0 : for(size_t i = 0; i < m_vNumberData.size(); ++i)
104 0 : base_type::update_subset_groups(m_vNumberData[i]);
105 0 : for(size_t i = 0; i < m_vBNDNumberData.size(); ++i)
106 0 : base_type::update_subset_groups(m_vBNDNumberData[i]);
107 0 : for(size_t i = 0; i < m_vVectorData.size(); ++i)
108 0 : base_type::update_subset_groups(m_vVectorData[i]);
109 0 : }
110 :
111 :
112 : ////////////////////////////////////////////////////////////////////////////////
113 : // assembling functions
114 : ////////////////////////////////////////////////////////////////////////////////
115 :
116 : template<typename TDomain>
117 : template<typename TElem, typename TFVGeom>
118 0 : void NeumannBoundaryFV<TDomain>::
119 : prep_elem_loop(const ReferenceObjectID roid, const int si)
120 : {
121 0 : update_subset_groups();
122 0 : m_si = si;
123 :
124 : // register subsetIndex at Geometry
125 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID,m_order);
126 :
127 : try{
128 0 : geo.update_local(roid, m_lfeID);
129 : }
130 0 : UG_CATCH_THROW("NeumannBoundaryFV::prep_elem_loop:"
131 : " Cannot update Finite Volume Geometry.");
132 :
133 : // request subset indices as boundary subset. This will force the
134 : // creation of boundary subsets when calling geo.update
135 :
136 0 : for(size_t i = 0; i < m_vNumberData.size(); ++i){
137 0 : if(!m_vNumberData[i].InnerSSGrp.contains(m_si)) continue;
138 0 : for(size_t s = 0; s < m_vNumberData[i].BndSSGrp.size(); ++s){
139 : const int si = m_vNumberData[i].BndSSGrp[s];
140 0 : geo.add_boundary_subset(si);
141 : }
142 : }
143 0 : for(size_t i = 0; i < m_vBNDNumberData.size(); ++i){
144 0 : if(!m_vBNDNumberData[i].InnerSSGrp.contains(m_si)) continue;
145 0 : for(size_t s = 0; s < m_vBNDNumberData[i].BndSSGrp.size(); ++s){
146 : const int si = m_vBNDNumberData[i].BndSSGrp[s];
147 0 : geo.add_boundary_subset(si);
148 : }
149 : }
150 0 : for(size_t i = 0; i < m_vVectorData.size(); ++i){
151 0 : if(!m_vVectorData[i].InnerSSGrp.contains(m_si)) continue;
152 0 : for(size_t s = 0; s < m_vVectorData[i].BndSSGrp.size(); ++s){
153 : const int si = m_vVectorData[i].BndSSGrp[s];
154 0 : geo.add_boundary_subset(si);
155 : }
156 : }
157 :
158 : // clear imports, since we will set them afterwards
159 : this->clear_imports();
160 :
161 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
162 :
163 : // set lin defect fct for imports
164 0 : for(size_t data = 0; data < m_vNumberData.size(); ++data)
165 : {
166 0 : if(!m_vNumberData[data].InnerSSGrp.contains(m_si)) continue;
167 0 : m_vNumberData[data].import.set_fct(id,
168 : &m_vNumberData[data],
169 : &NumberData::template lin_def<TElem, TFVGeom>);
170 :
171 0 : this->register_import(m_vNumberData[data].import);
172 : m_vNumberData[data].import.set_rhs_part();
173 : }
174 0 : }
175 :
176 : template<typename TDomain>
177 : template<typename TElem, typename TFVGeom>
178 0 : void NeumannBoundaryFV<TDomain>::
179 : prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[])
180 : {
181 : // update Geometry for this element
182 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID,m_order);
183 : try{
184 0 : geo.update(elem,vCornerCoords,&(this->subset_handler()));
185 : }
186 0 : UG_CATCH_THROW("NeumannBoundaryFV::prep_elem: "
187 : "Cannot update Finite Volume Geometry.");
188 :
189 0 : for(size_t i = 0; i < m_vNumberData.size(); ++i)
190 0 : if(m_vNumberData[i].InnerSSGrp.contains(m_si))
191 0 : m_vNumberData[i].template extract_bip<TElem, TFVGeom>(geo);
192 0 : }
193 :
194 : template<typename TDomain>
195 : template<typename TElem, typename TFVGeom>
196 0 : void NeumannBoundaryFV<TDomain>::
197 : add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[])
198 : {
199 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID,m_order);
200 : typedef typename TFVGeom::BF BF;
201 :
202 : // Number Data
203 0 : for(size_t data = 0; data < m_vNumberData.size(); ++data){
204 0 : if(!m_vNumberData[data].InnerSSGrp.contains(m_si)) continue;
205 : size_t ip = 0;
206 0 : for(size_t s = 0; s < m_vNumberData[data].BndSSGrp.size(); ++s){
207 : const int si = m_vNumberData[data].BndSSGrp[s];
208 : const std::vector<BF>& vBF = geo.bf(si);
209 :
210 0 : for(size_t b = 0; b < vBF.size(); ++b){
211 0 : const int co = vBF[b].node_id();
212 :
213 0 : for(size_t i = 0; i < vBF[b].num_ip(); ++i, ++ip){
214 0 : d(_C_, co) -= m_vNumberData[data].import[ip]
215 0 : * vBF[b].volume() * vBF[b].weight(i);
216 : }
217 : }
218 : }
219 : }
220 :
221 : // conditional Number Data
222 0 : for(size_t data = 0; data < m_vBNDNumberData.size(); ++data){
223 0 : if(!m_vBNDNumberData[data].InnerSSGrp.contains(m_si)) continue;
224 0 : for(size_t s = 0; s < m_vBNDNumberData[data].BndSSGrp.size(); ++s) {
225 : const int si = m_vBNDNumberData[data].BndSSGrp[s];
226 : const std::vector<BF>& vBF = geo.bf(si);
227 :
228 0 : for(size_t b = 0; b < vBF.size(); ++b){
229 0 : number val = 0.0;
230 0 : const int co = vBF[b].node_id();
231 :
232 0 : for(size_t i = 0; i < vBF[b].num_ip(); ++i){
233 0 : if(!(*m_vBNDNumberData[data].functor)(val, vBF[b].global_ip(i), this->time(), si))
234 0 : continue;
235 :
236 0 : d(_C_, co) -= val * vBF[b].volume() * vBF[b].weight(i);
237 : }
238 : }
239 : }
240 : }
241 :
242 : // vector data
243 0 : for(size_t data = 0; data < m_vVectorData.size(); ++data){
244 0 : if(!m_vVectorData[data].InnerSSGrp.contains(m_si)) continue;
245 0 : for(size_t s = 0; s < m_vVectorData[data].BndSSGrp.size(); ++s){
246 : const int si = m_vVectorData[data].BndSSGrp[s];
247 : const std::vector<BF>& vBF = geo.bf(si);
248 :
249 0 : for(size_t b = 0; b < vBF.size(); ++b){
250 : MathVector<dim> val;
251 0 : const int co = vBF[b].node_id();
252 :
253 0 : for(size_t i = 0; i < vBF[b].num_ip(); ++i){
254 0 : (*m_vVectorData[data].functor)(val, vBF[b].global_ip(i), this->time(), si);
255 :
256 0 : d(_C_, co) -= vBF[b].weight(i) * VecDot(val, vBF[b].normal());
257 : }
258 : }
259 : }
260 : }
261 0 : }
262 :
263 : template<typename TDomain>
264 : template<typename TElem, typename TFVGeom>
265 0 : void NeumannBoundaryFV<TDomain>::
266 : finish_elem_loop()
267 : {
268 : // remove subsetIndex from Geometry
269 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID,m_order);
270 :
271 : // unrequest subset indices as boundary subset. This will force the
272 : // creation of boundary subsets when calling geo.update
273 :
274 0 : for(size_t i = 0; i < m_vNumberData.size(); ++i){
275 0 : if(!m_vNumberData[i].InnerSSGrp.contains(m_si)) continue;
276 0 : for(size_t s = 0; s < m_vNumberData[i].BndSSGrp.size(); ++s){
277 : const int si = m_vNumberData[i].BndSSGrp[s];
278 0 : geo.remove_boundary_subset(si);
279 : geo.reset_curr_elem();
280 : }
281 : }
282 :
283 0 : for(size_t i = 0; i < m_vBNDNumberData.size(); ++i){
284 0 : if(!m_vBNDNumberData[i].InnerSSGrp.contains(m_si)) continue;
285 0 : for(size_t s = 0; s < m_vBNDNumberData[i].BndSSGrp.size(); ++s){
286 : const int si = m_vBNDNumberData[i].BndSSGrp[s];
287 0 : geo.remove_boundary_subset(si);
288 : geo.reset_curr_elem();
289 : }
290 : }
291 :
292 0 : for(size_t i = 0; i < m_vVectorData.size(); ++i){
293 0 : if(!m_vVectorData[i].InnerSSGrp.contains(m_si)) continue;
294 0 : for(size_t s = 0; s < m_vVectorData[i].BndSSGrp.size(); ++s){
295 : const int si = m_vVectorData[i].BndSSGrp[s];
296 0 : geo.remove_boundary_subset(si);
297 : geo.reset_curr_elem();
298 : }
299 : }
300 0 : }
301 :
302 : ////////////////////////////////////////////////////////////////////////////////
303 : // Number Data
304 : ////////////////////////////////////////////////////////////////////////////////
305 :
306 : template<typename TDomain>
307 : template<typename TElem, typename TFVGeom>
308 0 : void NeumannBoundaryFV<TDomain>::NumberData::
309 : lin_def(const LocalVector& u,
310 : std::vector<std::vector<number> > vvvLinDef[],
311 : const size_t nip)
312 : {
313 : // get finite volume geometry
314 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(This->m_lfeID,This->m_order);
315 : typedef typename TFVGeom::BF BF;
316 :
317 : size_t ip = 0;
318 0 : for(size_t s = 0; s < this->BndSSGrp.size(); ++s)
319 : {
320 : const int si = this->BndSSGrp[s];
321 : const std::vector<BF>& vBF = geo.bf(si);
322 0 : for(size_t i = 0; i < vBF.size(); ++i){
323 0 : const int co = vBF[i].node_id();
324 0 : vvvLinDef[ip][_C_][co] -= vBF[i].volume();
325 : }
326 : }
327 0 : }
328 :
329 : template<typename TDomain>
330 : template<typename TElem, typename TFVGeom>
331 0 : void NeumannBoundaryFV<TDomain>::NumberData::
332 : extract_bip(const TFVGeom& geo)
333 : {
334 : typedef typename TFVGeom::BF BF;
335 : static const int locDim = TElem::dim;
336 :
337 : std::vector<MathVector<locDim> >* vLocIP = local_ips<locDim>();
338 :
339 : vLocIP->clear();
340 : vGloIP.clear();
341 0 : for(size_t s = 0; s < this->BndSSGrp.size(); s++)
342 : {
343 : const int si = this->BndSSGrp[s];
344 : const std::vector<BF>& vBF = geo.bf(si);
345 0 : for(size_t i = 0; i < vBF.size(); ++i)
346 : {
347 : const BF& bf = vBF[i];
348 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip){
349 0 : vLocIP->push_back(bf.local_ip(ip));
350 0 : vGloIP.push_back(bf.global_ip(ip));
351 : }
352 : }
353 : }
354 :
355 0 : import.template set_local_ips<locDim>(&(*vLocIP)[0], vLocIP->size());
356 0 : import.set_global_ips(&vGloIP[0], vGloIP.size());
357 0 : }
358 :
359 : template<typename TDomain>
360 : template <int refDim>
361 : std::vector<MathVector<refDim> >* NeumannBoundaryFV<TDomain>::NumberData::local_ips()
362 : {
363 : switch (refDim)
364 : {
365 0 : case 1: return (std::vector<MathVector<refDim> >*)(&vLocIP_dim1);
366 0 : case 2: return (std::vector<MathVector<refDim> >*)(&vLocIP_dim2);
367 0 : case 3: return (std::vector<MathVector<refDim> >*)(&vLocIP_dim3);
368 : }
369 : }
370 :
371 : ////////////////////////////////////////////////////////////////////////////////
372 : // register assemble functions
373 : ////////////////////////////////////////////////////////////////////////////////
374 :
375 : #ifdef UG_DIM_1
376 : template<>
377 0 : void NeumannBoundaryFV<Domain1d>::register_all_funcs(int order)
378 : {
379 : // Edge
380 0 : register_func<RegularEdge, DimFVGeometry<dim, 1> >();
381 0 : }
382 : #endif
383 :
384 : #ifdef UG_DIM_2
385 : template<>
386 0 : void NeumannBoundaryFV<Domain2d>::register_all_funcs(int order)
387 : {
388 : // Edge
389 0 : register_func<RegularEdge, DimFVGeometry<dim, 1> >();
390 :
391 : // Triangle
392 0 : switch(order)
393 : {
394 0 : case 1: {typedef FVGeometry<1, Triangle, dim> FVGeom;
395 0 : register_func<Triangle, FVGeom >(); break;}
396 0 : case 2: {typedef FVGeometry<2, Triangle, dim> FVGeom;
397 0 : register_func<Triangle, FVGeom >(); break;}
398 0 : case 3: {typedef FVGeometry<3, Triangle, dim> FVGeom;
399 0 : register_func<Triangle, FVGeom >(); break;}
400 0 : default: {typedef DimFVGeometry<dim, 2> FVGeom;
401 0 : register_func<Triangle, FVGeom >(); break;}
402 : }
403 :
404 : // Quadrilateral
405 0 : switch(order) {
406 0 : case 1: {typedef FVGeometry<1, Quadrilateral, dim> FVGeom;
407 0 : register_func<Quadrilateral, FVGeom >(); break;}
408 0 : case 2: {typedef FVGeometry<2, Quadrilateral, dim> FVGeom;
409 0 : register_func<Quadrilateral, FVGeom >(); break;}
410 0 : case 3: {typedef FVGeometry<3, Quadrilateral, dim> FVGeom;
411 0 : register_func<Quadrilateral, FVGeom >(); break;}
412 0 : default: {typedef DimFVGeometry<dim, 2> FVGeom;
413 0 : register_func<Quadrilateral, FVGeom >(); break;}
414 : }
415 0 : }
416 : #endif
417 :
418 : #ifdef UG_DIM_3
419 : template<>
420 0 : void NeumannBoundaryFV<Domain3d>::register_all_funcs(int order)
421 : {
422 : // Edge
423 0 : register_func<RegularEdge, DimFVGeometry<dim, 1> >();
424 :
425 : // Tetrahedron
426 0 : switch(order)
427 : {
428 0 : case 1: {typedef FVGeometry<1, Tetrahedron, dim> FVGeom;
429 0 : register_func<Tetrahedron, FVGeom >(); break;}
430 0 : case 2: {typedef FVGeometry<2, Tetrahedron, dim> FVGeom;
431 0 : register_func<Tetrahedron, FVGeom >(); break;}
432 0 : case 3: {typedef FVGeometry<3, Tetrahedron, dim> FVGeom;
433 0 : register_func<Tetrahedron, FVGeom >(); break;}
434 0 : default: {typedef DimFVGeometry<dim, 3> FVGeom;
435 0 : register_func<Tetrahedron, FVGeom >(); break;}
436 : }
437 :
438 : // Prism
439 0 : switch(order) {
440 0 : case 1: {typedef FVGeometry<1, Prism, dim> FVGeom;
441 0 : register_func<Prism, FVGeom >(); break;}
442 0 : default: {typedef DimFVGeometry<dim, 3> FVGeom;
443 0 : register_func<Prism, FVGeom >(); break;}
444 : }
445 :
446 : // Hexahedron
447 0 : switch(order)
448 : {
449 0 : case 1: {typedef FVGeometry<1, Hexahedron, dim> FVGeom;
450 0 : register_func<Hexahedron, FVGeom >(); break;}
451 0 : case 2: {typedef FVGeometry<2, Hexahedron, dim> FVGeom;
452 0 : register_func<Hexahedron, FVGeom >(); break;}
453 0 : case 3: {typedef FVGeometry<3, Hexahedron, dim> FVGeom;
454 0 : register_func<Hexahedron, FVGeom >(); break;}
455 0 : default: {typedef DimFVGeometry<dim, 3> FVGeom;
456 0 : register_func<Hexahedron, FVGeom >(); break;}
457 : }
458 0 : }
459 : #endif
460 :
461 : template<typename TDomain>
462 : template<typename TElem, typename TFVGeom>
463 0 : void NeumannBoundaryFV<TDomain>::register_func()
464 : {
465 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
466 : typedef this_type T;
467 :
468 : this->clear_add_fct(id);
469 :
470 : this->set_prep_elem_loop_fct(id, &T::template prep_elem_loop<TElem, TFVGeom>);
471 : this->set_prep_elem_fct( id, &T::template prep_elem<TElem, TFVGeom>);
472 : this->set_add_rhs_elem_fct( id, &T::template add_rhs_elem<TElem, TFVGeom>);
473 : this->set_fsh_elem_loop_fct( id, &T::template finish_elem_loop<TElem, TFVGeom>);
474 :
475 : this->set_add_jac_A_elem_fct( id, &T::template add_jac_A_elem<TElem, TFVGeom>);
476 : this->set_add_jac_M_elem_fct( id, &T::template add_jac_M_elem<TElem, TFVGeom>);
477 : this->set_add_def_A_elem_fct( id, &T::template add_def_A_elem<TElem, TFVGeom>);
478 : this->set_add_def_M_elem_fct( id, &T::template add_def_M_elem<TElem, TFVGeom>);
479 0 : }
480 :
481 : ////////////////////////////////////////////////////////////////////////////////
482 : // explicit template instantiations
483 : ////////////////////////////////////////////////////////////////////////////////
484 :
485 : #ifdef UG_DIM_1
486 : template class NeumannBoundaryFV<Domain1d>;
487 : #endif
488 : #ifdef UG_DIM_2
489 : template class NeumannBoundaryFV<Domain2d>;
490 : #endif
491 : #ifdef UG_DIM_3
492 : template class NeumannBoundaryFV<Domain3d>;
493 : #endif
494 :
495 : } // namespace ug
496 :
|
