Loading [MathJax]/extensions/tex2jax.js
Plugins
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
filter_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2014: G-CSC, Goethe University Frankfurt
3 * Author: Christian Wehner
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#ifndef __H__UG__NAVIER_STOKES__INCOMPRESSIBLE__FILTER__IMPL__
34#define __H__UG__NAVIER_STOKES__INCOMPRESSIBLE__FILTER__IMPL__
35
36#include "common/util/provider.h"
37
38namespace ug{
39namespace NavierStokes{
40
41template<typename VType,typename TElem,typename TDomain,typename TGridFunction>
42void averageByVolume(PeriodicAttachmentAccessor<TElem,Attachment<VType> >& acUHat,
43 PeriodicAttachmentAccessor<TElem,Attachment<number> >& acVolume,TDomain& domain,SmartPtr<TGridFunction> uInfo,std::vector<WallObject<TGridFunction> > walls){
44 typedef typename TDomain::grid_type TGrid;
45 typedef typename TGridFunction::template traits<TElem>::const_iterator ElemIterator;
46 PeriodicBoundaryManager* pbm = (domain.grid())->periodic_boundary_manager();
47 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
48 {
49 for (int j=0;j<(int)walls.size();j++) if ((int)si == (int)walls[j].si()) continue;
50 ElemIterator elemIter = uInfo->template begin<TElem>(si);
51 ElemIterator elemIterEnd = uInfo->template end<TElem>(si);
52 for( ;elemIter !=elemIterEnd; elemIter++)
53 {
54 TElem* elem = *elemIter;
55 if (pbm && pbm->is_slave(elem)) continue;
56 acUHat[elem]/=(number)acVolume[elem];
57 }
58 }
59}
60
61template <typename vector_t>
62inline
63typename vector_t::value_type
64VecDistanceMaxNorm(const vector_t& v1, const vector_t& v2)
65{
66 typename vector_t::value_type m = std::abs(v1[0]-v2[0]);
67 typedef typename vector_t::size_type size_type;
68
69 for(size_type i = 1; i < v1.size(); ++i)
70 {
71 m = std::max(m, std::abs(v1[i]-v2[i]) );
72 }
73 return m;
74}
75
76// compute average element size of grid
77template<typename TGridFunction>
78number ConstantBoxFilter<TGridFunction>::compute_average_element_size(SmartPtr<TGridFunction> u){
79 typedef typename TGridFunction::const_element_iterator const_iterator;
80 typedef typename TGridFunction::element_type element_type;
81 typename TGridFunction::domain_type::position_accessor_type& aaPos
82 = u->domain()->position_accessor();
83 number averageElemSize = 0;
84 std::vector<MathVector<dim> > vCorner;
85 size_t nrOfElem = 0;
86 const_iterator iter = m_uInfo->template begin<elem_type>();
87 const_iterator iterEnd = m_uInfo->template end<elem_type>();
88 for(; iter != iterEnd; ++iter)
89 {
90 elem_type* elem = *iter;
91 ReferenceObjectID roid = elem->reference_object_id();
92 CollectCornerCoordinates(vCorner, *elem, aaPos);
93 averageElemSize += ElementSize<dim>(roid, &vCorner[0]);
94 nrOfElem++;
95 }
96 return averageElemSize/(number)nrOfElem;
97}
98
99template <typename TGridFunction>
100template <typename VType>
101void ConstantBoxFilter<TGridFunction>::collectSides(PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acUHat,
102 PeriodicAttachmentAccessor<side_type,Attachment<number> >& acVolume,
103 std::vector< MathVector<dim> >& coord,
104 VType values[DimFV1Geometry<dim>::maxNumSCV],
105 number volumes[DimFV1Geometry<dim>::maxNumSCV],
106 std::vector<side_type*>& sides
107 ){
108 domain_type& domain = *m_uInfo->domain().get();
109 PeriodicBoundaryManager* pbm = (domain.grid())->periodic_boundary_manager();
110 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
111 // get position accessor
112 typedef typename domain_type::position_accessor_type position_accessor_type;
113 const position_accessor_type& posAcc = domain.position_accessor();
114 side_secure_container assoSides;
115 size_t startIndex = sides.size()-1;
116 // handle start side
117 side_type* s = sides[startIndex];
118 MathVector<dim> bary = posAcc[s->vertex(0)];
119 for (size_t k=1;k<s->num_vertices();k++){
120 bary += posAcc[s->vertex(k)];
121 }
122 bary /= s->num_vertices();
123 number width=0.5*m_width;
124 for (size_t w=0;w<m_walls.size();w++){
125 width = std::min(width,m_walls[w].dist(bary));
126 }
127 for (size_t k=0;k<coord.size();k++){
128 number dist = VecDistanceMaxNorm(coord[k],bary);
129 if (dist<width){
130 if ((!pbm)||(pbm && !pbm->is_slave(s))){
131 acVolume[s]+=volumes[k];
132 acUHat[s]+=values[k];
133 }
134 }
135 }
136 for (size_t i=startIndex;i<sides.size();i++){
137 for (size_t v=0;v<sides[i]->num_vertices();v++){
138 domain.grid()->associated_elements(assoSides,sides[i]->vertex(v));
139 for (size_t j=0;j<assoSides.size();j++){
140 s = assoSides[j];
141 // check if new side
142 bool isNew = true;
143 for (int k=sides.size()-1;k>=0;k--){
144 if (sides[k]==s){
145 isNew=false;
146 break;
147 }
148 }
149 if (isNew==false) continue;
150 // compute dof coord
151 bary = posAcc[s->vertex(0)];
152 for (size_t k=1;k<s->num_vertices();k++){
153 bary += posAcc[s->vertex(k)];
154 }
155 bary /= s->num_vertices();
156 number width=0.5*m_width;
157 for (size_t w=0;w<m_walls.size();w++){
158 width = std::min(width,m_walls[w].dist(bary));
159 }
160 // check if inside filter region
161 bool inside=false;
162 for (size_t k=0;k<coord.size();k++){
163 number dist = VecDistanceMaxNorm(coord[k],bary);
164 if (dist<0.5*m_width) inside=true;
165 if (dist<width){
166 if ((!pbm)||(pbm && !pbm->is_slave(s))){
167 acVolume[s]+=volumes[k];
168 acUHat[s]+=values[k];
169 }
170 }
171
172 }
173 // if inside add to list
174 if (inside==true) sides.push_back(s);
175 }
176 }
177 }
178}
179
180// constant box filter for side data
181template <typename TGridFunction>
182template <typename VType>
183void ConstantBoxFilter<TGridFunction>::apply_filter(PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acUHat,
184 SmartPtr<TGridFunction> u,PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acU){
185 bool useGridFunction;
186 if (u==SPNULL) useGridFunction = false; else useGridFunction = true;
187 // set attachments to zero
188 domain_type& domain = *m_uInfo->domain().get();
189 SetAttachmentValues(m_acSideVolume, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
190 SetAttachmentValues(acUHat, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
191
192 std::vector<DoFIndex> multInd;
193
194 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
195
196 // get position accessor
197 typedef typename domain_type::position_accessor_type position_accessor_type;
198 const position_accessor_type& posAcc = domain.position_accessor();
199
200 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
201 side_secure_container elemsides;
202
203 DimFV1Geometry<dim> geo;
204 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
205 {
206 // get iterators
207 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
208 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
209
210 // loop elements of dimension
211 for( ;iter !=iterEnd; ++iter)
212 {
213 // get Elem
214 elem_type* elem = *iter;
215
216 // get corners of element
217 size_t noc=elem->num_vertices();
218 for(size_t i = 0; i < noc; ++i)
219 coCoord[i] = posAcc[elem->vertex(i)];
220
221 // get sides
222 domain.grid()->associated_elements_sorted(elemsides, static_cast<elem_type*>(elem) );
223
224 // evaluate finite volume geometry
225 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
226
227 // reference object id
228 ReferenceObjectID roid = elem->reference_object_id();
229
230 // get trial space
231 const LocalShapeFunctionSet<dim>& rTrialSpace =
232 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::CROUZEIX_RAVIART,dim, 1));
233
234 VType values[DimFV1Geometry<dim>::maxNumSCV];
235 number volumes[DimFV1Geometry<dim>::maxNumSCV];
236 std::vector< MathVector<dim> > baryCoords(geo.num_scv());
237 std::vector<side_type*> sides;
238 sides.reserve(20);
239 std::vector< Region<dim> > regions(1);
240 regions[0].firstIndex=0;
241 regions[0].periodicOffset=0;
242
243
244 for (size_t i=0;i<geo.num_scv();i++){
245 typename DimFV1Geometry<dim>::SCV scv = geo.scv(i);
246 baryCoords.resize(geo.num_scv());
247 MathVector<dim> scvLocalBary = scv.local_corner(0);
248 baryCoords[i] = scv.global_corner(0);
249 for (size_t co=1;co<scv.num_corners();co++){
250 scvLocalBary+=scv.local_corner(co);
251 baryCoords[i]+=scv.global_corner(co);
252 }
253 scvLocalBary/=scv.num_corners();
254 baryCoords[i]/=scv.num_corners();
255
256 std::vector<number> vShape;
257 rTrialSpace.shapes(vShape, scvLocalBary);
258 //MathVector<dim> localValue = 0;
259
260 size_t nos=elemsides.size();
261 // sum up shapes
262 values[i] = 0;
263 for (size_t sh=0;sh<nos;sh++){
264 VType localValue;
265 if (useGridFunction)
266 for (size_t d=0;d<dim;d++){
267 u->dof_indices(elemsides[sh], d, multInd);
268 assignVal(localValue,d,DoFRef(*u,multInd[0]));
269 }
270 else
271 localValue = acU[elemsides[sh]];
272 localValue *= vShape[sh];
273 values[i] += localValue;
274 }
275 values[i] *= scv.volume();
276 volumes[i] = scv.volume();
277 }
278 // check neighborhood for sides in the filter region
279 sides.push_back(elemsides[0]);
280 collectSides(acUHat,m_acSideVolume,baryCoords,values,volumes,sides);
281 PeriodicBoundaryManager* pbm = (domain.grid())->periodic_boundary_manager();
282 // check for periodic sides and add corresponding side
283 if (pbm){
284 for (size_t i=0;i<sides.size();i++){
285 side_type* s = sides[i];
286 if ( (pbm->is_slave(s)) || (pbm->is_master(s)) ) {
287 side_type* sNew;
288 if (pbm->is_slave(s))
289 sNew = pbm->master(s);
290 else
291 sNew = (*pbm->slaves(s))[0];
292 bool isNew = true;
293 for (size_t j=0;j<sides.size();j++){
294 if (sides[j]==sNew){
295 isNew=false;
296 break;
297 }
298 }
299 if (isNew==true){
300 sides.push_back(sNew);
301 // identify base region
302 size_t regInd;
303 if (regions.size()==0)
304 regInd=0;
305 else{
306 for (regInd=regions.size()-1; ;regInd--){
307 if (i>regions[regInd].firstIndex) break;
308 if (regInd==0) break;
309 }
310 }
311 // compute periodic offset
312 MathVector<dim> sReg0Co = posAcc[s->vertex(0)];
313 MathVector<dim> sReg1Co = posAcc[sNew->vertex(0)];
314 for (size_t k=1;k<sNew->num_vertices();k++){
315 sReg0Co += posAcc[s->vertex(k)];
316 sReg1Co += posAcc[sNew->vertex(k)];
317 }
318 sReg0Co/=s->num_vertices();
319 sReg1Co/=sNew->num_vertices();
320 MathVector<dim> pOffset;
321 for (size_t k=0;k<dim;k++){
322 pOffset[k] = regions[regInd].periodicOffset[k]+sReg0Co[k]-sReg1Co[k];
323 }
324 // add region
325 size_t rsize = regions.size();
326 regions.resize(rsize+1);
327 regions[rsize].firstIndex=sides.size()-1;
328 regions[rsize].periodicOffset=pOffset;
329 std::vector< MathVector<dim> > baryCoordsNew(geo.num_scv());
330 for (size_t k=0;k<baryCoordsNew.size();k++){
331 for (size_t d0=0;d0<dim;d0++){
332 baryCoordsNew[k][d0] = baryCoords[k][d0]-pOffset[d0];
333 }
334 }
335 // find associated sides
336 collectSides(acUHat,m_acSideVolume,baryCoordsNew,values,volumes,sides);
337 }
338 }
339 }
340 }
341 }
342 }
343 averageByVolume(acUHat,m_acSideVolume,domain,m_uInfo,m_walls);
344 if (useGridFunction) copyWallData(acUHat,u,m_walls);
345 else copyWallData(acUHat,acU,m_uInfo,m_walls);
346 if (dim==2) constrainingSideAveraging<TGridFunction,side_type,ConstrainingEdge,VType>(acUHat,m_uInfo);
347 else {
348 constrainingSideAveraging<TGridFunction,side_type,ConstrainingTriangle,VType>(acUHat,m_uInfo);
349 constrainingSideAveraging<TGridFunction,side_type,ConstrainingQuadrilateral,VType>(acUHat,m_uInfo);
350 }
351};
352
353// compute average element size of grid
354template<typename TGridFunction>
355number VariableBoxFilter<TGridFunction>::compute_average_element_size(SmartPtr<TGridFunction> u){
356 typedef typename TGridFunction::const_element_iterator const_iterator;
357 typedef typename TGridFunction::element_type element_type;
358 typename TGridFunction::domain_type::position_accessor_type& aaPos
359 = u->domain()->position_accessor();
360 number averageElemSize = 0;
361 std::vector<MathVector<dim> > vCorner;
362 size_t nrOfElem = 0;
363 const_iterator iter = m_uInfo->template begin<elem_type>();
364 const_iterator iterEnd = m_uInfo->template end<elem_type>();
365 for(; iter != iterEnd; ++iter)
366 {
367 elem_type* elem = *iter;
368 ReferenceObjectID roid = elem->reference_object_id();
369 CollectCornerCoordinates(vCorner, *elem, aaPos);
370 averageElemSize += ElementSize<dim>(roid, &vCorner[0]);
371 nrOfElem++;
372 }
373 return averageElemSize/(number)nrOfElem;
374}
375
376template <typename TGridFunction>
377template <typename VType>
378void VariableBoxFilter<TGridFunction>::collectSides(PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acUHat,
379 PeriodicAttachmentAccessor<side_type,Attachment<number> >& acVolume,
380 std::vector< MathVector<dim> >& coord,
381 VType values[DimFV1Geometry<dim>::maxNumSCV],
382 number volumes[DimFV1Geometry<dim>::maxNumSCV],
383 std::vector<side_type*>& sides
384 ){
385 domain_type& domain = *m_uInfo->domain().get();
386 PeriodicBoundaryManager* pbm = (domain.grid())->periodic_boundary_manager();
387 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
388 // get position accessor
389 std::vector<DoFIndex> multInd;
390 typedef typename domain_type::position_accessor_type position_accessor_type;
391 const position_accessor_type& posAcc = domain.position_accessor();
392 side_secure_container assoSides;
393 size_t startIndex = sides.size()-1;
394 // handle start side
395 side_type* s = sides[startIndex];
396 MathVector<dim> bary = posAcc[s->vertex(0)];
397 for (size_t k=1;k<s->num_vertices();k++){
398 bary += posAcc[s->vertex(k)];
399 }
400 bary /= s->num_vertices();
401 m_width->inner_dof_indices(s,0, multInd);
402 number localwidth=0.5*DoFRef(*m_width,multInd[0]);
403 for (size_t w=0;w<m_walls.size();w++){
404 localwidth = std::min(localwidth,m_walls[w].dist(bary));
405 }
406 for (size_t k=0;k<coord.size();k++){
407 number dist = VecDistanceMaxNorm(coord[k],bary);
408 if (dist<0.5*localwidth){
409 if ((!pbm)||(pbm && !pbm->is_slave(s))){
410 acVolume[s]+=volumes[k];
411 acUHat[s]+=values[k];
412 }
413 }
414 }
415 for (size_t i=startIndex;i<sides.size();i++){
416 for (size_t v=0;v<sides[i]->num_vertices();v++){
417 domain.grid()->associated_elements(assoSides,sides[i]->vertex(v));
418 for (size_t j=0;j<assoSides.size();j++){
419 s = assoSides[j];
420 // check if new side
421 bool isNew = true;
422 for (int k=sides.size()-1;k>=0;k--){
423 if (sides[k]==s){
424 isNew=false;
425 break;
426 }
427 }
428 if (isNew==false) continue;
429 // compute dof coord
430 bary = posAcc[s->vertex(0)];
431 for (size_t k=1;k<s->num_vertices();k++){
432 bary += posAcc[s->vertex(k)];
433 }
434 bary /= s->num_vertices();
435 // get local filterwidth
436 m_width->inner_dof_indices(s,0, multInd);
437 localwidth=0.5*DoFRef(*m_width,multInd[0]);
438 for (size_t w=0;w<m_walls.size();w++){
439 localwidth = std::min(localwidth,m_walls[w].dist(bary));
440 }
441 // check if inside filter region
442 bool inside=false;
443 for (size_t k=0;k<coord.size();k++){
444 number dist = VecDistanceMaxNorm(coord[k],bary);
445 if (dist<0.5*localwidth){
446 if ((!pbm)||(pbm && !pbm->is_slave(s))){
447 acVolume[s]+=volumes[k];
448 acUHat[s]+=values[k];
449 }
450 }
451 if (dist<0.5*m_maxwidth) inside=true;
452 }
453 // if inside add to list
454 if (inside==true) sides.push_back(s);
455 }
456 }
457 }
458}
459
460// variable box filter for side data
461template <typename TGridFunction>
462template <typename VType>
463void VariableBoxFilter<TGridFunction>::apply_filter(PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acUHat,
464 SmartPtr<TGridFunction> u,PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acU){
465 bool useGridFunction;
466 if (u==SPNULL) useGridFunction = false; else useGridFunction = true;
467 // set attachments to zero
468 domain_type& domain = *m_uInfo->domain().get();
469 SetAttachmentValues(m_acSideVolume, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
470 SetAttachmentValues(acUHat, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
471
472 std::vector<DoFIndex> multInd;
473
474 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
475
476 // get position accessor
477 typedef typename domain_type::position_accessor_type position_accessor_type;
478 const position_accessor_type& posAcc = domain.position_accessor();
479
480 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
481 side_secure_container elemsides;
482
483 DimFV1Geometry<dim> geo;
484 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
485 {
486 // get iterators
487 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
488 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
489
490 // loop elements of dimension
491 for( ;iter !=iterEnd; ++iter)
492 {
493 // get Elem
494 elem_type* elem = *iter;
495
496 // get corners of element
497 size_t noc=elem->num_vertices();
498 for(size_t i = 0; i < noc; ++i)
499 coCoord[i] = posAcc[elem->vertex(i)];
500
501 // get sides
502 domain.grid()->associated_elements_sorted(elemsides, static_cast<elem_type*>(elem) );
503
504 // evaluate finite volume geometry
505 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
506
507 // reference object id
508 ReferenceObjectID roid = elem->reference_object_id();
509
510 // get trial space
511 const LocalShapeFunctionSet<dim>& rTrialSpace =
512 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::CROUZEIX_RAVIART,dim, 1));
513
514 VType values[DimFV1Geometry<dim>::maxNumSCV];
515 number volumes[DimFV1Geometry<dim>::maxNumSCV];
516 std::vector< MathVector<dim> > baryCoords(geo.num_scv());
517 std::vector<side_type*> sides;
518 sides.reserve(20);
519 std::vector< Region<dim> > regions(1);
520 regions[0].firstIndex=0;
521 regions[0].periodicOffset=0;
522
523
524 for (size_t i=0;i<geo.num_scv();i++){
525 typename DimFV1Geometry<dim>::SCV scv = geo.scv(i);
526 baryCoords.resize(geo.num_scv());
527 MathVector<dim> scvLocalBary = scv.local_corner(0);
528 baryCoords[i] = scv.global_corner(0);
529 for (size_t co=1;co<scv.num_corners();co++){
530 scvLocalBary+=scv.local_corner(co);
531 baryCoords[i]+=scv.global_corner(co);
532 }
533 scvLocalBary/=scv.num_corners();
534 baryCoords[i]/=scv.num_corners();
535
536 std::vector<number> vShape;
537 rTrialSpace.shapes(vShape, scvLocalBary);
538 //MathVector<dim> localValue = 0;
539
540 size_t nos=elemsides.size();
541 // sum up shapes
542 values[i] = 0;
543 for (size_t sh=0;sh<nos;sh++){
544 VType localValue;
545 if (useGridFunction)
546 for (size_t d=0;d<dim;d++){
547 u->dof_indices(elemsides[sh], d, multInd);
548 assignVal(localValue,d,DoFRef(*u,multInd[0]));
549 }
550 else
551 localValue = acU[elemsides[sh]];
552 localValue *= vShape[sh];
553 values[i] += localValue;
554 }
555 values[i] *= scv.volume();
556 volumes[i] = scv.volume();
557 }
558 // check neighborhood for sides in the filter region
559 sides.push_back(elemsides[0]);
560 collectSides(acUHat,m_acSideVolume,baryCoords,values,volumes,sides);
561 PeriodicBoundaryManager* pbm = (domain.grid())->periodic_boundary_manager();
562 // check for periodic sides and add corresponding side
563 if (pbm){
564 for (size_t i=0;i<sides.size();i++){
565 side_type* s = sides[i];
566 if ( (pbm->is_slave(s)) || (pbm->is_master(s)) ) {
567 side_type* sNew;
568 if (pbm->is_slave(s))
569 sNew = pbm->master(s);
570 else
571 sNew = (*pbm->slaves(s))[0];
572 bool isNew = true;
573 for (size_t j=0;j<sides.size();j++){
574 if (sides[j]==sNew){
575 isNew=false;
576 break;
577 }
578 }
579 if (isNew==true){
580 sides.push_back(sNew);
581 // identify base region
582 size_t regInd;
583 if (regions.size()==0)
584 regInd=0;
585 else{
586 for (regInd=regions.size()-1; ;regInd--){
587 if (i>regions[regInd].firstIndex) break;
588 if (regInd==0) break;
589 }
590 }
591 // compute periodic offset
592 MathVector<dim> sReg0Co = posAcc[s->vertex(0)];
593 MathVector<dim> sReg1Co = posAcc[sNew->vertex(0)];
594 for (size_t k=1;k<sNew->num_vertices();k++){
595 sReg0Co += posAcc[s->vertex(k)];
596 sReg1Co += posAcc[sNew->vertex(k)];
597 }
598 sReg0Co/=s->num_vertices();
599 sReg1Co/=sNew->num_vertices();
600 MathVector<dim> pOffset;
601 for (size_t k=0;k<dim;k++){
602 pOffset[k] = regions[regInd].periodicOffset[k]+sReg0Co[k]-sReg1Co[k];
603 }
604 // add region
605 size_t rsize = regions.size();
606 regions.resize(rsize+1);
607 regions[rsize].firstIndex=sides.size()-1;
608 regions[rsize].periodicOffset=pOffset;
609 std::vector< MathVector<dim> > baryCoordsNew(geo.num_scv());
610 for (size_t k=0;k<baryCoordsNew.size();k++){
611 for (size_t d0=0;d0<dim;d0++){
612 baryCoordsNew[k][d0] = baryCoords[k][d0]-pOffset[d0];
613 }
614 }
615 // find associated sides
616 collectSides(acUHat,m_acSideVolume,baryCoordsNew,values,volumes,sides);
617 }
618 }
619 }
620 }
621 }
622 }
623 averageByVolume(acUHat,m_acSideVolume,domain,m_uInfo,m_walls);
624 if (useGridFunction) copyWallData(acUHat,u,m_walls);
625 else copyWallData(acUHat,acU,m_uInfo,m_walls);
626 if (dim==2) constrainingSideAveraging<TGridFunction,side_type,ConstrainingEdge,VType>(acUHat,m_uInfo);
627 else {
628 constrainingSideAveraging<TGridFunction,side_type,ConstrainingTriangle,VType>(acUHat,m_uInfo);
629 constrainingSideAveraging<TGridFunction,side_type,ConstrainingQuadrilateral,VType>(acUHat,m_uInfo);
630 }
631};
632
633// fv1 Filter for vertex data
634template <typename TGridFunction>
635template <typename VType>
636void FV1BoxFilter<TGridFunction>::apply_filter(PeriodicAttachmentAccessor<vertex_type,Attachment<VType> >& acUHat,
637 SmartPtr<TGridFunction> u,PeriodicAttachmentAccessor<vertex_type,Attachment<VType> >& acU){
638 bool useGridFunction = true;
639 if (u==SPNULL) useGridFunction = false;
640
641 // set attachments to zero
642 domain_type& domain = *m_uInfo->domain().get();
643 SetAttachmentValues(m_acVertexVolume, domain.grid()->template begin<vertex_type>(), domain.grid()->template end<vertex_type>(), 0);
644 SetAttachmentValues(acUHat, domain.grid()->template begin<vertex_type>(), domain.grid()->template end<vertex_type>(), 0);
645
646 std::vector<DoFIndex> multInd;
647
648 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
649
650 // get position accessor
651 typedef typename domain_type::position_accessor_type position_accessor_type;
652 const position_accessor_type& posAcc = domain.position_accessor();
653
654 DimFV1Geometry<dim> geo;
655 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
656 {
657 // get iterators
658 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
659 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
660
661 // loop elements of dimension
662 for( ;iter !=iterEnd; ++iter)
663 {
664 // get Elem
665 elem_type* elem = *iter;
666
667 // get corners of element
668 size_t noc=elem->num_vertices();
669 for(size_t i = 0; i < noc; ++i)
670 coCoord[i] = posAcc[elem->vertex(i)];
671
672 // evaluate finite volume geometry
673 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
674
675 // reference object id
676 ReferenceObjectID roid = elem->reference_object_id();
677
678 // get trial space
679 const LocalShapeFunctionSet<dim>& rTrialSpace =
680 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE,dim, 1));
681
682 for (size_t i=0;i<geo.num_scv();i++){
683 typename DimFV1Geometry<dim>::SCV scv = geo.scv(i);
684 MathVector<dim> scvLocalBary = scv.local_corner(0);
685 for (size_t co=1;co<scv.num_corners();co++){
686 scvLocalBary+=scv.local_corner(co);
687 }
688 scvLocalBary/=scv.num_corners();
689 size_t node = scv.node_id();
690
691 std::vector<number> vShape;
692 rTrialSpace.shapes(vShape, scvLocalBary);
693 // sum up shapes
694 VType value;
695 value = 0;
696 for (size_t sh=0;sh<noc;sh++){
697 VType localValue;
698 if (useGridFunction)
699 for (size_t d=0;d<dim;d++){
700 u->dof_indices(elem->vertex(sh), d, multInd);
701 assignVal(localValue,d,DoFRef(*u,multInd[0]));
702 }
703 else
704 localValue = acU[elem->vertex(sh)];
705 localValue *= vShape[sh];
706 value += localValue;
707 }
708 value *= scv.volume();
709 m_acVertexVolume[elem->vertex(node)] += scv.volume();
710 acUHat[elem->vertex(node)] += value;
711 }
712 }
713 }
714 averageByVolume(acUHat,m_acVertexVolume,domain,m_uInfo,m_walls);
715 if (useGridFunction) copyWallData(acUHat,u,m_walls);
716 else copyWallData(acUHat,acU,m_uInfo,m_walls);
717}
718
719// fv1 Filter for vertex data
720template <typename TGridFunction>
721void FV1BoxFilter<TGridFunction>::compute_filterwidth_fv1()
722{
723 // set attachments to zero
724 domain_type& domain = *m_uInfo->domain().get();
725 SetAttachmentValues(m_acVertexVolume, domain.grid()->template begin<vertex_type>(), domain.grid()->template end<vertex_type>(), 0);
726 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
727
728 // get position accessor
729 typedef typename domain_type::position_accessor_type position_accessor_type;
730 const position_accessor_type& posAcc = domain.position_accessor();
731
732 DimFV1Geometry<dim> geo;
733 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
734 {
735 // get iterators
736 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
737 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
738
739 // loop elements of dimension
740 for( ;iter !=iterEnd; ++iter)
741 {
742 // get Elem
743 elem_type* elem = *iter;
744
745 // get corners of element
746 size_t noc=elem->num_vertices();
747 for(size_t i = 0; i < noc; ++i)
748 coCoord[i] = posAcc[elem->vertex(i)];
749
750 // evaluate finite volume geometry
751 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
752
753 for (size_t i=0;i<geo.num_scv();i++){
754 typename DimFV1Geometry<dim>::SCV scv = geo.scv(i);
755 size_t node = scv.node_id();
756 m_acVertexVolume[elem->vertex(node)] += scv.volume();
757 }
758 }
759 }
760}
761
762// fv1 Filter for side data
763template <typename TGridFunction>
764void FV1BoxFilter<TGridFunction>::compute_filterwidth_fvcr()
765{
766 domain_type& domain = *m_uInfo->domain().get();
767 SetAttachmentValues(m_acSideVolume, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
768 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
769
770 // get position accessor
771 typedef typename domain_type::position_accessor_type position_accessor_type;
772 const position_accessor_type& posAcc = domain.position_accessor();
773
774 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
775 side_secure_container sides;
776
777 DimFV1Geometry<dim> geo;
778 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
779 {
780 // get iterators
781 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
782 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
783
784 // loop elements of dimension
785 for( ;iter !=iterEnd; ++iter)
786 {
787 // get Elem
788 elem_type* elem = *iter;
789
790 // get corners of element
791 size_t noc=elem->num_vertices();
792 for(size_t i = 0; i < noc; ++i)
793 coCoord[i] = posAcc[elem->vertex(i)];
794
795 // get sides
796 domain.grid()->associated_elements_sorted(sides, static_cast<elem_type*>(elem) );
797
798 // evaluate finite volume geometry
799 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
800
801 for (size_t i=0;i<geo.num_scv();i++){
802 typename DimFV1Geometry<dim>::SCV scv = geo.scv(i);
803 size_t node = scv.node_id();
804 // find sides associated with scv
805 for (size_t j=0;j<sides.size();j++){
806 side_type* sidej = sides[j];
807 for (size_t sco=0;sco<sidej->num_vertices();sco++){
808 if (VecDistance(coCoord[node],posAcc[sidej->vertex(sco)])<1e-12){
809 m_acSideVolume[sidej] += scv.volume();
810 }
811 }
812 }
813 }
814 }
815 }
816}
817
818
819// fv1 Filter for side data
820template <typename TGridFunction>
821template <typename VType>
822void FV1BoxFilter<TGridFunction>::apply_filter(PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acUHat,
823 SmartPtr<TGridFunction> u,PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acU){
824 bool useGridFunction = true;
825 if (u==SPNULL) useGridFunction = false;
826
827 // set attachments to zero
828 domain_type& domain = *m_uInfo->domain().get();
829 SetAttachmentValues(m_acSideVolume, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
830 SetAttachmentValues(acUHat, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
831
832 std::vector<DoFIndex> multInd;
833
834 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
835
836 // get position accessor
837 typedef typename domain_type::position_accessor_type position_accessor_type;
838 const position_accessor_type& posAcc = domain.position_accessor();
839
840 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
841 side_secure_container sides;
842
843 DimFV1Geometry<dim> geo;
844 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
845 {
846 // get iterators
847 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
848 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
849
850 // loop elements of dimension
851 for( ;iter !=iterEnd; ++iter)
852 {
853 // get Elem
854 elem_type* elem = *iter;
855
856 // get corners of element
857 size_t noc=elem->num_vertices();
858 for(size_t i = 0; i < noc; ++i)
859 coCoord[i] = posAcc[elem->vertex(i)];
860
861 // get sides
862 domain.grid()->associated_elements_sorted(sides, static_cast<elem_type*>(elem) );
863
864 // evaluate finite volume geometry
865 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
866
867 // reference object id
868 ReferenceObjectID roid = elem->reference_object_id();
869
870 // get trial space
871 const LocalShapeFunctionSet<dim>& rTrialSpace =
872 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::CROUZEIX_RAVIART,dim, 1));
873
874 for (size_t i=0;i<geo.num_scv();i++){
875 typename DimFV1Geometry<dim>::SCV scv = geo.scv(i);
876 MathVector<dim> scvLocalBary = scv.local_corner(0);
877 for (size_t co=1;co<scv.num_corners();co++){
878 scvLocalBary+=scv.local_corner(co);
879 }
880 scvLocalBary/=scv.num_corners();
881 size_t node = scv.node_id();
882
883 std::vector<number> vShape;
884 rTrialSpace.shapes(vShape, scvLocalBary);
885// MathVector<dim> localValue = 0;
886
887 size_t nos=sides.size();
888 // sum up shapes
889 VType value;
890 value = 0;
891 for (size_t sh=0;sh<nos;sh++){
892 VType localValue;
893 if (useGridFunction)
894 for (size_t d=0;d<dim;d++){
895 u->dof_indices(sides[sh], d, multInd);
896 assignVal(localValue,d,DoFRef(*u,multInd[0]));
897 }
898 else
899 localValue = acU[sides[sh]];
900 localValue *= vShape[sh];
901 value += localValue;
902 }
903 value *= scv.volume();
904 // find sides associated with scv
905 for (size_t j=0;j<nos;j++){
906 side_type* sidej = sides[j];
907 for (size_t sco=0;sco<sidej->num_vertices();sco++){
908 if (VecDistance(coCoord[node],posAcc[sidej->vertex(sco)])<1e-12){
909 m_acSideVolume[sidej] += scv.volume();
910 acUHat[sidej] += value;
911 }
912 }
913 }
914 }
915 }
916 }
917 averageByVolume(acUHat,m_acSideVolume,domain,m_uInfo,m_walls);
918 if (useGridFunction) copyWallData(acUHat,u,m_walls);
919 else copyWallData(acUHat,acU,m_uInfo,m_walls);
920 if (dim==2) constrainingSideAveraging<TGridFunction,side_type,ConstrainingEdge,VType>(acUHat,m_uInfo);
921 else {
922 constrainingSideAveraging<TGridFunction,side_type,ConstrainingTriangle,VType>(acUHat,m_uInfo);
923 constrainingSideAveraging<TGridFunction,side_type,ConstrainingQuadrilateral,VType>(acUHat,m_uInfo);
924 }
925}
926
927// fvcr Filter for side data
928template <typename TGridFunction>
929template <typename VType>
930void FVCRBoxFilter<TGridFunction>::apply_filter(PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acUHat,
931 SmartPtr<TGridFunction> u,PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acU){
932 bool useGridFunction = true;
933 if (u==SPNULL) useGridFunction = false;
934
935 // set attachments to zero
936 domain_type& domain = *m_uInfo->domain().get();
937 SetAttachmentValues(m_acSideVolume, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
938 SetAttachmentValues(acUHat, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
939
940 std::vector<DoFIndex> multInd;
941
942 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
943
944 // get position accessor
945 typedef typename domain_type::position_accessor_type position_accessor_type;
946 const position_accessor_type& posAcc = domain.position_accessor();
947
948 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
949 side_secure_container sides;
950
951 DimCRFVGeometry<dim> geo;
952 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
953 {
954 // get iterators
955 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
956 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
957
958 // loop elements of dimension
959 for( ;iter !=iterEnd; ++iter)
960 {
961 // get Elem
962 elem_type* elem = *iter;
963
964 // get corners of element
965 size_t noc=elem->num_vertices();
966 for(size_t i = 0; i < noc; ++i)
967 coCoord[i] = posAcc[elem->vertex(i)];
968
969 // get sides
970 domain.grid()->associated_elements_sorted(sides, static_cast<elem_type*>(elem) );
971
972 // evaluate finite volume geometry
973 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
974
975 // reference object id
976 ReferenceObjectID roid = elem->reference_object_id();
977
978 // get trial space
979 const LocalShapeFunctionSet<dim>& rTrialSpace =
980 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::CROUZEIX_RAVIART,dim, 1));
981
982 for (size_t i=0;i<geo.num_scv();i++){
983 typename DimCRFVGeometry<dim>::SCV scv = geo.scv(i);
984 MathVector<dim> scvLocalBary = scv.local_corner(0);
985 for (size_t co=1;co<scv.num_corners();co++){
986 scvLocalBary+=scv.local_corner(co);
987 }
988 scvLocalBary/=scv.num_corners();
989 size_t sidei = scv.node_id();
990
991 std::vector<number> vShape;
992 rTrialSpace.shapes(vShape, scvLocalBary);
993
994 size_t nos=sides.size();
995 // sum up shapes
996 VType value;
997 value = 0;
998 for (size_t sh=0;sh<nos;sh++){
999 VType localValue;
1000 if (useGridFunction)
1001 for (size_t d=0;d<dim;d++){
1002 u->dof_indices(sides[sh], d, multInd);
1003 assignVal(localValue,d,DoFRef(*u,multInd[0]));
1004 }
1005 else
1006 localValue = acU[sides[sh]];
1007 localValue *= vShape[sh];
1008 value += localValue;
1009 }
1010 value *= scv.volume();
1011 m_acSideVolume[sides[sidei]] += scv.volume();
1012 acUHat[sides[sidei]] += value;
1013 }
1014 }
1015 }
1016 averageByVolume(acUHat,m_acSideVolume,domain,m_uInfo,m_walls);
1017 if (useGridFunction) copyWallData(acUHat,u,m_walls);
1018 else copyWallData(acUHat,acU,m_uInfo,m_walls);
1019 if (dim==2) constrainingSideAveraging<TGridFunction,side_type,ConstrainingEdge,VType>(acUHat,m_uInfo);
1020 else {
1021 constrainingSideAveraging<TGridFunction,side_type,ConstrainingTriangle,VType>(acUHat,m_uInfo);
1022 constrainingSideAveraging<TGridFunction,side_type,ConstrainingQuadrilateral,VType>(acUHat,m_uInfo);
1023 }
1024}
1025
1026template <typename TGridFunction>
1027void FVCRBoxFilter<TGridFunction>::compute_filterwidth_fvcr()
1028{
1029 // set attachments to zero
1030 domain_type& domain = *m_uInfo->domain().get();
1031 SetAttachmentValues(m_acSideVolume, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
1032
1033 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
1034
1035 // get position accessor
1036 typedef typename domain_type::position_accessor_type position_accessor_type;
1037 const position_accessor_type& posAcc = domain.position_accessor();
1038
1039 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
1040 side_secure_container sides;
1041
1042 DimCRFVGeometry<dim> geo;
1043 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
1044 {
1045 // get iterators
1046 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
1047 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
1048
1049 // loop elements of dimension
1050 for( ;iter !=iterEnd; ++iter)
1051 {
1052 // get Elem
1053 elem_type* elem = *iter;
1054
1055 // get corners of element
1056 size_t noc=elem->num_vertices();
1057 for(size_t i = 0; i < noc; ++i)
1058 coCoord[i] = posAcc[elem->vertex(i)];
1059
1060 // get sides
1061 domain.grid()->associated_elements_sorted(sides, static_cast<elem_type*>(elem) );
1062
1063 // evaluate finite volume geometry
1064 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
1065
1066 for (size_t i=0;i<geo.num_scv();i++){
1067 typename DimCRFVGeometry<dim>::SCV scv = geo.scv(i);
1068 size_t sidei = scv.node_id();
1069 m_acSideVolume[sides[sidei]] += scv.volume();
1070 }
1071 }
1072 }
1073}
1074
1075template <typename TGridFunction>
1076void ElementBoxFilter<TGridFunction>::compute_filterwidth_fv1(){
1077 // set attachments to zero
1078 domain_type& domain = *m_uInfo->domain().get();
1079 SetAttachmentValues(m_acVertexVolume, domain.grid()->template begin<vertex_type>(), domain.grid()->template end<vertex_type>(), 0);
1080 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
1081
1082 // get position accessor
1083 typedef typename domain_type::position_accessor_type position_accessor_type;
1084 const position_accessor_type& posAcc = domain.position_accessor();
1085
1086 DimFV1Geometry<dim> geo;
1087 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
1088 {
1089 // get iterators
1090 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
1091 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
1092
1093 // loop elements of dimension
1094 for( ;iter !=iterEnd; ++iter)
1095 {
1096 // get Elem
1097 elem_type* elem = *iter;
1098
1099 // get corners of element
1100 size_t noc=elem->num_vertices();
1101 for(size_t i = 0; i < noc; ++i)
1102 coCoord[i] = posAcc[elem->vertex(i)];
1103
1104 // evaluate finite volume geometry
1105 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
1106
1107 for (size_t i=0;i<geo.num_scv();i++){
1108 typename DimFV1Geometry<dim>::SCV scv = geo.scv(i);
1109 for (size_t node=0;node<elem->num_vertices();node++){
1110 m_acVertexVolume[elem->vertex(node)] += scv.volume();
1111 }
1112 }
1113 }
1114 }
1115}
1116
1117// Element filter for vertex data
1118template <typename TGridFunction>
1119template <typename VType>
1120void ElementBoxFilter<TGridFunction>::apply_filter(PeriodicAttachmentAccessor<vertex_type,Attachment<VType> >& acUHat,
1121 SmartPtr<TGridFunction> u,PeriodicAttachmentAccessor<vertex_type,Attachment<VType> >& acU){
1122 bool useGridFunction = true;
1123 if (u==SPNULL) useGridFunction = false;
1124
1125 // set attachments to zero
1126 domain_type& domain = *m_uInfo->domain().get();
1127 SetAttachmentValues(m_acVertexVolume, domain.grid()->template begin<vertex_type>(), domain.grid()->template end<vertex_type>(), 0);
1128 SetAttachmentValues(acUHat, domain.grid()->template begin<vertex_type>(), domain.grid()->template end<vertex_type>(), 0);
1129
1130 std::vector<DoFIndex> multInd;
1131
1132 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
1133
1134 // get position accessor
1135 typedef typename domain_type::position_accessor_type position_accessor_type;
1136 const position_accessor_type& posAcc = domain.position_accessor();
1137
1138 DimFV1Geometry<dim> geo;
1139 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
1140 {
1141 // get iterators
1142 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
1143 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
1144
1145 // loop elements of dimension
1146 for( ;iter !=iterEnd; ++iter)
1147 {
1148 // get Elem
1149 elem_type* elem = *iter;
1150
1151 // get corners of element
1152 size_t noc=elem->num_vertices();
1153 for(size_t i = 0; i < noc; ++i)
1154 coCoord[i] = posAcc[elem->vertex(i)];
1155
1156 // evaluate finite volume geometry
1157 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
1158
1159 // reference object id
1160 ReferenceObjectID roid = elem->reference_object_id();
1161
1162 // get trial space
1163 const LocalShapeFunctionSet<dim>& rTrialSpace =
1164 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE,dim, 1));
1165
1166 for (size_t i=0;i<geo.num_scv();i++){
1167 typename DimFV1Geometry<dim>::SCV scv = geo.scv(i);
1168 MathVector<dim> scvLocalBary = scv.local_corner(0);
1169 for (size_t co=1;co<scv.num_corners();co++){
1170 scvLocalBary+=scv.local_corner(co);
1171 }
1172 scvLocalBary/=scv.num_corners();
1173
1174 std::vector<number> vShape;
1175 rTrialSpace.shapes(vShape, scvLocalBary);
1176 // sum up shapes
1177 VType value;
1178 value = 0;
1179 for (size_t sh=0;sh<noc;sh++){
1180 VType localValue;
1181 if (useGridFunction)
1182 for (size_t d=0;d<dim;d++){
1183 u->dof_indices(elem->vertex(sh), d, multInd);
1184 assignVal(localValue,d,DoFRef(*u,multInd[0]));
1185 }
1186 else
1187 localValue = acU[elem->vertex(sh)];
1188 localValue *= vShape[sh];
1189 value += localValue;
1190 }
1191 value *= scv.volume();
1192 for (size_t node=0;node<elem->num_vertices();node++){
1193 m_acVertexVolume[elem->vertex(node)] += scv.volume();
1194 acUHat[elem->vertex(node)] += value;
1195 }
1196 }
1197 }
1198 }
1199 averageByVolume(acUHat,m_acVertexVolume,domain,m_uInfo,m_walls);
1200 if (useGridFunction) copyWallData(acUHat,u,m_walls);
1201 else copyWallData(acUHat,acU,m_uInfo,m_walls);
1202}
1203
1204// fvcr Filter for side data
1205template <typename TGridFunction>
1206void ElementBoxFilter<TGridFunction>::compute_filterwidth_fvcr(){
1207 // set attachments to zero
1208 domain_type& domain = *m_uInfo->domain().get();
1209 SetAttachmentValues(m_acSideVolume, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
1210 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
1211
1212 // get position accessor
1213 typedef typename domain_type::position_accessor_type position_accessor_type;
1214 const position_accessor_type& posAcc = domain.position_accessor();
1215
1216 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
1217 side_secure_container sides;
1218
1219 DimCRFVGeometry<dim> geo;
1220 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
1221 {
1222 // get iterators
1223 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
1224 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
1225
1226 // loop elements of dimension
1227 for( ;iter !=iterEnd; ++iter)
1228 {
1229 // get Elem
1230 elem_type* elem = *iter;
1231
1232 // get corners of element
1233 size_t noc=elem->num_vertices();
1234 for(size_t i = 0; i < noc; ++i)
1235 coCoord[i] = posAcc[elem->vertex(i)];
1236
1237 // get sides
1238 domain.grid()->associated_elements(sides, static_cast<elem_type*>(elem) );
1239
1240 // evaluate finite volume geometry
1241 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
1242
1243 for (size_t i=0;i<geo.num_scv();i++){
1244 typename DimCRFVGeometry<dim>::SCV scv = geo.scv(i);
1245 // sum up shapes
1246 for (size_t j=0;j<sides.size();j++){
1247 m_acSideVolume[sides[j]] += scv.volume();
1248 }
1249 }
1250 }
1251 }
1252}
1253
1254
1255// fvcr Filter for side data
1256template <typename TGridFunction>
1257template <typename VType>
1258void ElementBoxFilter<TGridFunction>::apply_filter(PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acUHat,
1259 SmartPtr<TGridFunction> u,PeriodicAttachmentAccessor<side_type,Attachment<VType> >& acU){
1260 bool useGridFunction = true;
1261 if (u==SPNULL) useGridFunction = false;
1262
1263 // set attachments to zero
1264 domain_type& domain = *m_uInfo->domain().get();
1265 SetAttachmentValues(m_acSideVolume, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
1266 SetAttachmentValues(acUHat, domain.grid()->template begin<side_type>(), domain.grid()->template end<side_type>(), 0);
1267
1268 std::vector<DoFIndex> multInd;
1269
1270 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
1271
1272 // get position accessor
1273 typedef typename domain_type::position_accessor_type position_accessor_type;
1274 const position_accessor_type& posAcc = domain.position_accessor();
1275
1276 typedef typename grid_type::template traits<side_type>::secure_container side_secure_container;
1277 side_secure_container sides;
1278
1279 DimCRFVGeometry<dim> geo;
1280 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
1281 {
1282 // get iterators
1283 ElemIterator iter = m_uInfo->template begin<elem_type>(si);
1284 ElemIterator iterEnd = m_uInfo->template end<elem_type>(si);
1285
1286 // loop elements of dimension
1287 for( ;iter !=iterEnd; ++iter)
1288 {
1289 // get Elem
1290 elem_type* elem = *iter;
1291
1292 // get corners of element
1293 size_t noc=elem->num_vertices();
1294 for(size_t i = 0; i < noc; ++i)
1295 coCoord[i] = posAcc[elem->vertex(i)];
1296
1297 // get sides
1298 domain.grid()->associated_elements_sorted(sides, static_cast<elem_type*>(elem) );
1299
1300 // evaluate finite volume geometry
1301 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
1302
1303 // reference object id
1304 ReferenceObjectID roid = elem->reference_object_id();
1305
1306 // get trial space
1307 const LocalShapeFunctionSet<dim>& rTrialSpace =
1308 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::CROUZEIX_RAVIART,dim, 1));
1309
1310 for (size_t i=0;i<geo.num_scv();i++){
1311 typename DimCRFVGeometry<dim>::SCV scv = geo.scv(i);
1312 MathVector<dim> scvLocalBary = scv.local_corner(0);
1313 for (size_t co=1;co<scv.num_corners();co++){
1314 scvLocalBary+=scv.local_corner(co);
1315 }
1316 scvLocalBary/=scv.num_corners();
1317
1318 std::vector<number> vShape;
1319 rTrialSpace.shapes(vShape, scvLocalBary);
1320
1321 size_t nos=sides.size();
1322 // sum up shapes
1323 VType value;
1324 value = 0;
1325 for (size_t sh=0;sh<nos;sh++){
1326 VType localValue;
1327 if (useGridFunction)
1328 for (size_t d=0;d<dim;d++){
1329 u->dof_indices(sides[sh], d, multInd);
1330 assignVal(localValue,d,DoFRef(*u,multInd[0]));
1331 }
1332 else
1333 localValue = acU[sides[sh]];
1334 localValue *= vShape[sh];
1335 value += localValue;
1336 }
1337 value *= scv.volume();
1338 for (size_t j=0;j<sides.size();j++){
1339 m_acSideVolume[sides[j]] += scv.volume();
1340 acUHat[sides[j]] += value;
1341 }
1342 }
1343 }
1344 }
1345 averageByVolume(acUHat,m_acSideVolume,domain,m_uInfo,m_walls);
1346 if (useGridFunction) copyWallData(acUHat,u,m_walls);
1347 else copyWallData(acUHat,acU,m_uInfo,m_walls);
1348 if (dim==2) constrainingSideAveraging<TGridFunction,side_type,ConstrainingEdge,VType>(acUHat,m_uInfo);
1349 else {
1350 constrainingSideAveraging<TGridFunction,side_type,ConstrainingTriangle,VType>(acUHat,m_uInfo);
1351 constrainingSideAveraging<TGridFunction,side_type,ConstrainingQuadrilateral,VType>(acUHat,m_uInfo);
1352 }
1353}
1354
1355
1356
1357} // namespace NavierStokes
1358} // end namespace ug
1359
1360#endif /* __H__UG__NAVIER_STOKES__INCOMPRESSIBLE__FILTER__IMPL__ */
s
parameterString s
Definition Biogas.lua:2
NavierStokes
function NavierStokes(fcts, subsets, discType)
ug::DimCRFVGeometry::SCV::num_corners
size_t num_corners() const
ug::DimCRFVGeometry::SCV::node_id
size_t node_id() const
ug::DimCRFVGeometry::SCV::local_corner
const MathVector< dim > & local_corner(size_t co) const
ug::DimCRFVGeometry::num_scv
size_t num_scv() const
ug::DimCRFVGeometry::scv
const SCV & scv(size_t i) const
ug::DimCRFVGeometry::update
void update(GridObject *elem, const MathVector< worldDim > *vCornerCoords, const ISubsetHandler *ish=NULL)
ug::DimFV1Geometry::SCV::local_corner
const MathVector< dim > & local_corner(size_t co) const
ug::DimFV1Geometry::SCV::volume
number volume() const
ug::DimFV1Geometry::SCV::global_corner
const MathVector< worldDim > & global_corner(size_t co) const
ug::DimFV1Geometry::SCV::node_id
size_t node_id() const
ug::DimFV1Geometry::SCV::num_corners
size_t num_corners() const
ug::DimFV1Geometry::num_scv
size_t num_scv() const
ug::DimFV1Geometry::update
void update(GridObject *elem, const MathVector< worldDim > *vCornerCoords, const ISubsetHandler *ish=NULL)
ug::DimFV1Geometry::scv
const SCV & scv(size_t i) const
ug::LocalShapeFunctionSet::shapes
virtual void shapes(std::vector< std::vector< shape_type > > &vvShape, const std::vector< MathVector< dim > > &vLocPos) const=0
ug::MathVector::size
std::size_t size() const
ug::NavierStokes::ConstantBoxFilter::elem_type
TGridFunction::template dim_traits< dim >::grid_base_object elem_type
Definition filter.h:393
ug::NavierStokes::ConstantBoxFilter::position_accessor_type
domain_type::position_accessor_type position_accessor_type
Definition filter.h:408
ug::NavierStokes::ConstantBoxFilter::compute_average_element_size
number compute_average_element_size(SmartPtr< TGridFunction > u)
Definition filter_impl.h:78
ug::NavierStokes::ConstantBoxFilter::apply_filter
void apply_filter(PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaUHat, SmartPtr< TGridFunction > u, PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaU)
Definition filter_impl.h:183
ug::NavierStokes::ConstantBoxFilter::ElemIterator
TGridFunction::template dim_traits< dim >::const_iterator ElemIterator
element iterator
Definition filter.h:411
ug::NavierStokes::ConstantBoxFilter::collectSides
void collectSides(PeriodicAttachmentAccessor< side_type, Attachment< VType > > &acUHat, PeriodicAttachmentAccessor< side_type, Attachment< number > > &acVolume, std::vector< MathVector< dim > > &coord, VType values[DimFV1Geometry< dim >::maxNumSCV], number volumes[DimFV1Geometry< dim >::maxNumSCV], std::vector< side_type * > &sides)
Definition filter_impl.h:101
ug::NavierStokes::ConstantBoxFilter::side_type
elem_type::side side_type
Definition filter.h:399
ug::NavierStokes::ConstantBoxFilter::domain_type
TGridFunction::domain_type domain_type
Definition filter.h:396
ug::NavierStokes::ElementBoxFilter::elem_type
TGridFunction::template dim_traits< dim >::grid_base_object elem_type
Definition filter.h:958
ug::NavierStokes::ElementBoxFilter::compute_filterwidth_fvcr
void compute_filterwidth_fvcr()
Definition filter_impl.h:1206
ug::NavierStokes::ElementBoxFilter::apply_filter
void apply_filter(PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaUHat, SmartPtr< TGridFunction > u, PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaU)
Definition filter_impl.h:1258
ug::NavierStokes::ElementBoxFilter::side_type
elem_type::side side_type
Definition filter.h:964
ug::NavierStokes::ElementBoxFilter::compute_filterwidth_fv1
void compute_filterwidth_fv1()
Definition filter_impl.h:1076
ug::NavierStokes::ElementBoxFilter::domain_type
TGridFunction::domain_type domain_type
Definition filter.h:961
ug::NavierStokes::ElementBoxFilter::ElemIterator
TGridFunction::template dim_traits< dim >::const_iterator ElemIterator
Definition filter.h:976
ug::NavierStokes::ElementBoxFilter::position_accessor_type
domain_type::position_accessor_type position_accessor_type
Definition filter.h:973
ug::NavierStokes::FV1BoxFilter::side_type
elem_type::side side_type
Definition filter.h:729
ug::NavierStokes::FV1BoxFilter::compute_filterwidth_fvcr
void compute_filterwidth_fvcr()
Definition filter_impl.h:764
ug::NavierStokes::FV1BoxFilter::position_accessor_type
domain_type::position_accessor_type position_accessor_type
Definition filter.h:738
ug::NavierStokes::FV1BoxFilter::apply_filter
void apply_filter(PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaUHat, SmartPtr< TGridFunction > u, PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaU)
Definition filter_impl.h:822
ug::NavierStokes::FV1BoxFilter::domain_type
TGridFunction::domain_type domain_type
Definition filter.h:726
ug::NavierStokes::FV1BoxFilter::elem_type
TGridFunction::template dim_traits< dim >::grid_base_object elem_type
Definition filter.h:723
ug::NavierStokes::FV1BoxFilter::ElemIterator
TGridFunction::template dim_traits< dim >::const_iterator ElemIterator
Definition filter.h:741
ug::NavierStokes::FV1BoxFilter::compute_filterwidth_fv1
void compute_filterwidth_fv1()
Definition filter_impl.h:721
ug::NavierStokes::FVCRBoxFilter::elem_type
TGridFunction::template dim_traits< dim >::grid_base_object elem_type
Definition filter.h:841
ug::NavierStokes::FVCRBoxFilter::position_accessor_type
domain_type::position_accessor_type position_accessor_type
Definition filter.h:856
ug::NavierStokes::FVCRBoxFilter::apply_filter
void apply_filter(PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaUHat, SmartPtr< TGridFunction > u, PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaU)
Definition filter_impl.h:930
ug::NavierStokes::FVCRBoxFilter::side_type
elem_type::side side_type
Definition filter.h:847
ug::NavierStokes::FVCRBoxFilter::ElemIterator
TGridFunction::template dim_traits< dim >::const_iterator ElemIterator
Definition filter.h:859
ug::NavierStokes::FVCRBoxFilter::domain_type
TGridFunction::domain_type domain_type
Definition filter.h:844
ug::NavierStokes::FVCRBoxFilter::compute_filterwidth_fvcr
void compute_filterwidth_fvcr()
Definition filter_impl.h:1027
ug::NavierStokes::VariableBoxFilter::compute_average_element_size
number compute_average_element_size(SmartPtr< TGridFunction > u)
Definition filter_impl.h:355
ug::NavierStokes::VariableBoxFilter::domain_type
TGridFunction::domain_type domain_type
Definition filter.h:531
ug::NavierStokes::VariableBoxFilter::ElemIterator
TGridFunction::template dim_traits< dim >::const_iterator ElemIterator
element iterator
Definition filter.h:546
ug::NavierStokes::VariableBoxFilter::collectSides
void collectSides(PeriodicAttachmentAccessor< side_type, Attachment< VType > > &acUHat, PeriodicAttachmentAccessor< side_type, Attachment< number > > &acVolume, std::vector< MathVector< dim > > &coord, VType values[DimFV1Geometry< dim >::maxNumSCV], number volumes[DimFV1Geometry< dim >::maxNumSCV], std::vector< side_type * > &sides)
Definition filter_impl.h:378
ug::NavierStokes::VariableBoxFilter::side_type
elem_type::side side_type
Definition filter.h:534
ug::NavierStokes::VariableBoxFilter::elem_type
TGridFunction::template dim_traits< dim >::grid_base_object elem_type
Definition filter.h:528
ug::NavierStokes::VariableBoxFilter::position_accessor_type
domain_type::position_accessor_type position_accessor_type
Definition filter.h:543
ug::NavierStokes::VariableBoxFilter::apply_filter
void apply_filter(PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaUHat, SmartPtr< TGridFunction > u, PeriodicAttachmentAccessor< side_type, Attachment< VType > > &aaU)
Definition filter_impl.h:463
ug::NavierStokes::WallObject
Definition filter.h:51
ug::PeriodicBoundaryManager::is_master
bool is_master(TElem *) const
ug::PeriodicBoundaryManager::slaves
Group< TElem >::SlaveContainer * slaves(TElem *e) const
ug::PeriodicBoundaryManager::master
TElem * master(TElem *e) const
ug::PeriodicBoundaryManager::is_slave
bool is_slave(TElem *) const
element_type
int element_type() const
ug::CollectCornerCoordinates
void CollectCornerCoordinates(int base_object_id, std::vector< typename TDomain::position_type > &vCornerCoordsOut, GridObject &elem, const TDomain &domain, bool clearContainer)
num_subsets
size_t num_subsets() const
ug::SetAttachmentValues
void SetAttachmentValues(TAttachmentAccessor &aaVal, TIter elemsBegin, TIter elemsEnd, const TVal &val)
SPNULL
const NullSmartPtr SPNULL
number
double number
ug::VecDistance
vector_t::value_type VecDistance(const vector_t &v1, const vector_t &v2)
ug::NavierStokes::VecDistanceMaxNorm
vector_t::value_type VecDistanceMaxNorm(const vector_t &v1, const vector_t &v2)
Definition filter_impl.h:64
ug::NavierStokes::averageByVolume
void averageByVolume(PeriodicAttachmentAccessor< TElem, Attachment< VType > > &acUHat, PeriodicAttachmentAccessor< TElem, Attachment< number > > &acVolume, TDomain &domain, SmartPtr< TGridFunction > uInfo, std::vector< WallObject< TGridFunction > > walls)
Definition filter_impl.h:42
ug::ReferenceObjectID
ReferenceObjectID
ug::DoFRef
const number & DoFRef(const TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
ug::domain_traits