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