Plugins
Loading...
Searching...
No Matches
pressure_separation.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2015: 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__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__FVCR__PRESSURE_SEPARATION__
34#define __H__UG__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__FVCR__PRESSURE_SEPARATION__
35
36#include "common/common.h"
37
49
50#ifdef UG_FOR_LUA
52#endif
53
54namespace ug{
55namespace NavierStokes{
56
57template<typename TElem,typename TPos,int dim>
58void computeElemBarycenter(MathVector<dim>& bary,TElem* elem,TPos posA){
59 bary = 0;
60 size_t noc=elem->num_vertices();
61 for (size_t i=0;i<noc;i++){
62 bary+=posA[elem->vertex(i)];
63 }
64 bary/=noc;
65}
66
67bool multMatVec(const std::vector<number>& avec, const std::vector<number>& b,std::vector<number>& c,size_t m,size_t n){
68 size_t i;
69 size_t j;
70 size_t count=0;
71 for(i = 0; i < m; i = i + 1){
72 c[i]=0.0;
73 for(j = 0; j < n; j = j + 1){
74 c[i] = c[i] + avec[count] * b[j];
75 count++;
76 }
77 }
78 return true;
79}
80
81
82bool solveLS(std::vector<number>& x,/* solution */
83 const std::vector<number>& matField,/* matrix given as field */
84 const std::vector<number>& b /* rhs */){
85 size_t n=b.size();
86 std::vector<std::vector<number> > P(n);
87 for(size_t i = 0; i < n; ++i) P[i].resize(n);
88
89 std::vector<number> Pvec(n*n);
90 std::vector<number> b2(n);
91
92 std::vector<std::vector<number> > A(n);
93 for(size_t i = 0; i < n; ++i) A[i].resize(n);
94 std::vector<std::vector<number> > L(n);
95 for(size_t i = 0; i < n; ++i) L[i].resize(n);
96 std::vector<std::vector<number> > U(n);
97 for(size_t i = 0; i < n; ++i) U[i].resize(n);
98// number A[n][n];
99// number L[n][n];
100// number U[n][n];
101
102 std::vector<number> z(n);
103 size_t i,k,j,l;
104 bool boolean=false;
105 number d,max;
106 size_t count=0;
107 for (i=0;i<n;i++){
108 for (j=0;j<n;j++){
109 A[i][j]=matField[count];
110 if (i!=j)
111 P[i][j]=0;
112 else
113 P[i][j]=1;
114 L[i][j]=0;
115 U[i][j]=0;
116 count ++;
117 }
118 }
119 for (i=0;i<n;i++){
120 j=i;
121 // Suche groesstes Element in Spalten
122 max=std::abs(A[i][i]);
123 for (k=i+1;k<n;k++){
124 if (std::abs(A[k][i])>max){
125 j=k;
126 max=std::abs(A[k][i]);
127 };
128 };
129 //debug
131 if (max<1e-14){
132 // A nicht invertierbar
133 // printf(".\n");
134 return false;
135 };
136 if (i!=j){
137 // Vertauschen der Zeilen
138 if (boolean==false){
139 // beim ersten Durchlauf Vertauschen der Spalten der
140 // Permutationsmatrix und der Zeilen von A,L.
141 boolean=true;
142 for (k=0;k<n;k++){
143 d=P[k][i];
144 P[k][i]=P[k][j];
145 P[k][j]=d;
146 d=A[i][k];
147 A[i][k]=A[j][k];
148 A[j][k]=d;
149 d=L[i][k];
150 L[i][k]=L[j][k];
151 L[j][k]=d;
152 };
153 } else{
154 for (k=0;k<n;k++){
155 // Vertauschen der Zeilen von P,A,L
156 d=P[i][k];
157 P[i][k]=P[j][k];
158 P[j][k]=d;
159 d=A[i][k];
160 A[i][k]=A[j][k];
161 A[j][k]=d;
162 d=L[i][k];
163 L[i][k]=L[j][k];
164 L[j][k]=d;
165 };
166 };
167 };
168 // Elimination
169 L[i][i]=1;
170 for (k=i+1;k<n;k++){
171 L[k][i]=A[k][i]/A[i][i];
172 for (l=i+1;l<n;l++){
173 A[k][l]=A[k][l]-L[k][i]*A[i][l];
174 };
175 };
176 };
177 // U ist der obere Dreiecksteil von A
178 for (i=0;i<n;i++){
179 for (j=i;j<n;j++){
180 U[i][j]=A[i][j];
181 };
182 };
183 // A*x = b <-> P*A*x = P*b
184 count = 0;
185 for (i=0;i<n;i++){
186 for (j=0;j<n;j++){
187 Pvec[count]=P[i][j];
188 count++;
189 }
190 }
191 multMatVec(Pvec,b,b2,n,n);
192 // Loese zuerst L*z = P*b
193 for (i=0;i<n;i++){
194 number s=b2[i];
195 for (j=0;j<i;j++){
196 s=s-z[j]*L[i][j];
197 };
198 z[i]=s;
199 };
200 // Loese U*x = z
201 for (i=n-1; ;i--){
202 number s=z[i];
203 for (j=n-1;j>i;j--){
204 s=s-x[j]*U[i][j];
205 };
206 x[i]=s/U[i][i];
207 if (i==0) break;
208 };
209 return true;
210};
211
212
213bool leastSquares(std::vector<number>& x,const std::vector<number>& mField,const std::vector<number>& b){
214 size_t m = b.size();
215 size_t n = x.size();
216 if (mField.size()!=m*n){
217 n = mField.size()/b.size();
218 }
219 std::vector<number> tmmField(n*n);
220 std::vector<number> tmb(n);
221 number z;
222 // compute A^t * A
223 for (size_t i=0;i<n;i++){
224 for (size_t j=i;j<n;j++){
225 z=0;
226 for (size_t k=0;k<m;k++){
227 z+=mField[n*k+i]*mField[n*k+j];
228 }
229 tmmField[j*n+i]=z;
230 if (j!=i) tmmField[i*n+j]=z;
231 }
232 tmb[i]=0;
233 for (size_t k=0;k<m;k++) tmb[i]+= mField[n*k+i]*b[k];
234 }
235 return solveLS(x,tmmField,tmb);
236}
237
238int bubblesort(std::vector<int>& list,std::vector<number> array){
239 int i,j,mini,k;
240 int length=array.size();
241 number min;
242 for (i=0;i<length;i++)
243 list[i]=i;
244 for (i=0;i<length-1;i++){
245 min = array[i];
246 mini= i;
247 for (j=i+1;j<length;j++){
248 if (array[j]<min){
249 mini = j;
250 min = array[j];
251 };
252 };
253 array[mini] = array[i];
254 k=list[i];
255 list[i] = list[mini];
256 list[mini] = k;
257 array[i] = min;
258 };
259 return 0;
260};
261
265template <typename TGridFunction>
267: public StdUserData<SeparatedPressureSource<TGridFunction>, MathVector<TGridFunction::dim>, TGridFunction::dim>,
268 virtual public INewtonUpdate
269 {
271 typedef typename TGridFunction::domain_type domain_type;
272
274 typedef typename TGridFunction::algebra_type algebra_type;
275
277 typedef typename domain_type::position_accessor_type position_accessor_type;
278
280 static const int dim = domain_type::dim;
281
283 typedef typename domain_type::grid_type grid_type;
284
286 typedef typename TGridFunction::template dim_traits<dim>::grid_base_object elem_type;
287
289 // typedef MathVector<dim> vecDim;
290 // typedef Attachment<vecDim> AMathVectorDim;
291
295
297 typedef typename TGridFunction::template dim_traits<dim>::const_iterator ElemIterator;
298
300 typedef typename TGridFunction::template traits<Vertex>::const_iterator VertexIterator;
301
302 private:
303 // old pressure attachment accessor
306
307 // pressure attachment accessor (interpolated pressure in vertices)
310
311 // volume attachment accessor
314
315 // level set grid function
317
318 // approximation space for level and surface grid
320
321 // grid
323
324 // pressure index
325 static const size_t _P_ = dim;
326
327 private:
328
331
332 public:
334
336 {
337 m_imSource = data;
338 }
339
341 {
343 for (int i=0;i<dim;i++){
344 f->set_entry(i, f_x);
345 }
346 set_source(f);
347 }
348
349 void set_source(number f_x, number f_y)
350 {
351 if (dim!=2){
352 UG_THROW("NavierStokes: Setting source vector of dimension 2"
353 " to a Discretization for world dim " << dim);
354 } else {
356 f->set_entry(0, f_x);
357 f->set_entry(1, f_y);
358 set_source(f);
359 }
360 }
361
362 void set_source(number f_x, number f_y, number f_z)
363 {
364 if (dim<3){
365 UG_THROW("NavierStokes: Setting source vector of dimension 3"
366 " to a Discretization for world dim " << dim);
367 }
368 else
369 {
371 f->set_entry(0, f_x);
372 f->set_entry(1, f_y);
373 f->set_entry(2, f_z);
374 set_source(f);
375 }
376 }
377
378#ifdef UG_FOR_LUA
379 void set_source(const char* fctName)
380 {
382 }
383#endif
384
385 public:
388 m_u = spGridFct;
389 domain_type& domain = *m_u->domain().get();
390 grid_type& grid = *domain.grid();
391 m_grid = &grid;
392 m_spApproxSpace = approxSpace;
393 set_source(0.0);
394 grid.template attach_to<elem_type>(m_aPOld);
395 grid.template attach_to<Vertex>(m_aP);
396 grid.template attach_to<Vertex>(m_aVol);
400 // set all values to zero
401 SetAttachmentValues(m_vol, m_u->template begin<Vertex>(), m_u->template end<Vertex>(), 0);
402 SetAttachmentValues(m_p, m_u->template begin<Vertex>(), m_u->template end<Vertex>(), 0);
403 SetAttachmentValues(m_pOld, m_u->template begin<elem_type>(), m_u->template end<elem_type>(), 0);
404 this->update();
405 }
406
408
409 template <int refDim>
410 inline void evaluate(MathVector<dim> vValue[],
411 const MathVector<dim> vGlobIP[],
412 number time, int si,
413 GridObject* elem,
414 const MathVector<dim> vCornerCoords[],
415 const MathVector<refDim> vLocIP[],
416 const size_t nip,
417 LocalVector* u,
418 const MathMatrix<refDim, dim>* vJT = NULL) const
419 {
420 UG_ASSERT(dynamic_cast<elem_type*>(elem) != NULL, "Unsupported element type");
421 elem_type* element = static_cast<elem_type*>(elem);
422
423 // reference object id
425
426 const size_t numVertices = element->num_vertices();
427 // get domain of grid function
428 const domain_type& domain = *m_u->domain().get();
429
430 // get position accessor
431 typedef typename domain_type::position_accessor_type position_accessor_type;
432 const position_accessor_type& posAcc = domain.position_accessor();
433
434// position_accessor_type aaPos = m_u->domain()->position_accessor();
435
436 // coord and vertex array
439 DimCRFVGeometry<dim> crfvgeo;
440
442
443 for(size_t i = 0; i < numVertices; ++i){
444 vVrt[i] = element->vertex(i);
445 coCoord[i] = posAcc[vVrt[i]];
446 };
447
448 for (size_t i=0;i<domain_traits<dim>::MaxNumVerticesOfElem;i++)
449 grad[i]*=0.0;
450
451 // evaluate finite volume geometry
452 crfvgeo.update(elem, &(coCoord[0]), domain.subset_handler().get());
453
454 // Lagrange 1 trial space
455 const LocalShapeFunctionSet<dim>& lagrange1 =
456 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
457
458 std::vector<number> shapes;
459
460 for (size_t ip=0;ip<crfvgeo.num_scvf();ip++){
461 const typename DimCRFVGeometry<dim>::SCVF& scvf = crfvgeo.scvf(ip);
462 number pinter=0;
463 lagrange1.shapes(shapes,scvf.local_ip());
464 for (size_t sh=0;sh<numVertices;sh++)
465 pinter += m_p[vVrt[sh]]*shapes[sh];
466 for (int d=0;d<dim;d++){
467 number flux = pinter*scvf.normal()[d];
468 grad[scvf.from()][d]-=flux;
469 grad[scvf.to()][d]+=flux;
470 }
471 }
472
473 for (size_t i=0;i<crfvgeo.num_scv();i++){
474 grad[i]/=crfvgeo.scv(i).volume();
475 }
476
477 // evaluate source data
478 (*m_imSource)(vValue,
479 vGlobIP,
480 time, si,
481 elem,
482 vCornerCoords,
483 vLocIP,
484 nip,
485 u,
486 vJT);
487 for (size_t i=0;i<nip;i++)
488 vValue[i]+=grad[i];
489 }; // evaluate
490
491 void update(){
492 // get domain
493 domain_type& domain = *m_u->domain().get();
494 // create Multiindex
495 std::vector<DoFIndex> multInd;
497
498 // coord and vertex array
501
502 // get position accessor
503 typedef typename domain_type::position_accessor_type position_accessor_type;
504 const position_accessor_type& posAcc = domain.position_accessor();
505
506 // set volume and p values to zero
507 SetAttachmentValues(m_vol, m_u->template begin<Vertex>(), m_u->template end<Vertex>(), 0);
508 SetAttachmentValues(m_p, m_u->template begin<Vertex>(), m_u->template end<Vertex>(), 0);
509
510 // set p^{old} = p^{old} + p^{sep}
511 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si){
512 ElemIterator iter = m_u->template begin<elem_type>(si);
513 ElemIterator iterEnd = m_u->template end<elem_type>(si);
514 for( ;iter !=iterEnd; ++iter)
515 {
516 elem_type* elem = *iter;
517 m_u->inner_dof_indices(elem, _P_, multInd);
518 m_pOld[elem]+=DoFRef(*m_u,multInd[0]);
519 }
520 }
521 // compute pressure in vertices by averaging
522 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si){
523 ElemIterator iter = m_u->template begin<elem_type>(si);
524 ElemIterator iterEnd = m_u->template end<elem_type>(si);
525 for( ;iter !=iterEnd; ++iter)
526 {
527 elem_type* elem = *iter;
528 number pValue = m_pOld[elem];
529 const size_t numVertices = elem->num_vertices();
530 for(size_t i = 0; i < numVertices; ++i){
531 vVrt[i] = elem->vertex(i);
532 coCoord[i] = posAcc[vVrt[i]];
533 };
534 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
535 for(size_t i = 0; i < numVertices; ++i){
536 number scvVol = geo.scv(i).volume();
537 m_vol[vVrt[i]]+=scvVol;
538 m_p[vVrt[i]]+=scvVol*pValue;
539 }
540 }
541 }
542 PeriodicBoundaryManager* pbm = (domain.grid())->periodic_boundary_manager();
543 // go over all vertices and average
544 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si){
545 VertexIterator iter = m_u->template begin<Vertex>(si);
546 VertexIterator iterEnd = m_u->template end<Vertex>(si);
547 for( ;iter !=iterEnd; ++iter)
548 {
549 Vertex* vrt = *iter;
550 if (pbm && pbm->is_slave(vrt)) continue;
551 m_p[vrt]/=m_vol[vrt];
552 }
553 }
554 }
555
556 private:
557 static const size_t max_number_of_ips = 20;
558
559 public:
561 const MathVector<dim>& globIP,
562 number time, int si) const
563 {
564 UG_THROW("LevelSetUserData: Need element.");
565 }
566
567 virtual void operator() (MathVector<dim> vValue[],
568 const MathVector<dim> vGlobIP[],
569 number time, int si, const size_t nip) const
570 {
571 UG_THROW("LevelSetUserData: Need element.");
572 }
573
574 virtual void compute(LocalVector* u, GridObject* elem,
575 const MathVector<dim> vCornerCoords[], bool bDeriv = false)
576 {
577 const int si = this->subset();
578 for(size_t s = 0; s < this->num_series(); ++s)
579 evaluate<dim>(this->values(s), this->ips(s), this->time(s), si,
580 elem, NULL, this->template local_ips<dim>(s),
581 this->num_ip(s), u);
582 }
583
585 const MathVector<dim> vCornerCoords[], bool bDeriv = false)
586 {
587 const int si = this->subset();
588 for(size_t s = 0; s < this->num_series(); ++s)
589 evaluate<dim>(this->values(s), this->ips(s), this->time(s), si,
590 elem, NULL, this->template local_ips<dim>(s),
591 this->num_ip(s), &(u->solution(this->time_point(s))));
592 }
593
595 virtual bool continuous() const {return false;}
596
598 virtual bool requires_grid_fct() const {return true;}
599};
600
601template <typename TGridFunction>
604 virtual public INewtonUpdate
605 {
607 typedef typename TGridFunction::domain_type domain_type;
608
610 typedef typename TGridFunction::algebra_type algebra_type;
611
613 typedef typename domain_type::position_accessor_type position_accessor_type;
614
616 static const int dim = domain_type::dim;
617
619 typedef typename domain_type::grid_type grid_type;
620
622 typedef typename TGridFunction::template dim_traits<dim>::grid_base_object elem_type;
623
625 // typedef MathVector<dim> vecDim;
626 // typedef Attachment<vecDim> AMathVectorDim;
627
631
633 typedef typename TGridFunction::template dim_traits<dim>::const_iterator ElemIterator;
634
636 typedef typename TGridFunction::template traits<Vertex>::const_iterator VertexIterator;
637
638 private:
639 // old pressure attachment accessor
642
643 // pressure attachment accessor (interpolated pressure in vertices)
646
647 // volume attachment accessor
650
651 // level set grid function
653
654 // approximation space for level and surface grid
656
657 // grid
659
660 // pressure index
661 static const size_t _P_ = dim;
662
663 private:
664
667
668 public:
670
672 {
673 m_imSource = data;
674 }
675
677 {
679 for (int i=0;i<dim;i++){
680 f->set_entry(i, f_x);
681 }
682 set_source(f);
683 }
684
685 void set_source(number f_x, number f_y)
686 {
687 if (dim!=2){
688 UG_THROW("NavierStokes: Setting source vector of dimension 2"
689 " to a Discretization for world dim " << dim);
690 } else {
692 f->set_entry(0, f_x);
693 f->set_entry(1, f_y);
694 set_source(f);
695 }
696 }
697
698 void set_source(number f_x, number f_y, number f_z)
699 {
700 if (dim<3){
701 UG_THROW("NavierStokes: Setting source vector of dimension 3"
702 " to a Discretization for world dim " << dim);
703 }
704 else
705 {
707 f->set_entry(0, f_x);
708 f->set_entry(1, f_y);
709 f->set_entry(2, f_z);
710 set_source(f);
711 }
712 }
713
714#ifdef UG_FOR_LUA
715 void set_source(const char* fctName)
716 {
717 set_source(LuaUserDataFactory<MathVector<dim>, dim>::create(fctName));
718 }
719#endif
720
721 public:
724 m_u = spGridFct;
725 domain_type& domain = *m_u->domain().get();
726 grid_type& grid = *domain.grid();
727 m_grid = &grid;
728 m_spApproxSpace = approxSpace;
729 set_source(0.0);
730 grid.template attach_to<elem_type>(m_aPOld);
731 grid.template attach_to<Vertex>(m_aP);
732 grid.template attach_to<Vertex>(m_aVol);
733 m_pOld.access(grid,m_aPOld);
734 m_p.access(grid,m_aP);
735 m_vol.access(grid,m_aVol);
736 // set all values to zero
737 SetAttachmentValues(m_vol, m_u->template begin<Vertex>(), m_u->template end<Vertex>(), 0);
738 SetAttachmentValues(m_p, m_u->template begin<Vertex>(), m_u->template end<Vertex>(), 0);
739 SetAttachmentValues(m_pOld, m_u->template begin<elem_type>(), m_u->template end<elem_type>(), 0);
740 this->update();
741 }
742
744
745 template <int refDim>
746 inline void evaluate(MathVector<dim> vValue[],
747 const MathVector<dim> vGlobIP[],
748 number time, int si,
749 GridObject* elem,
750 const MathVector<dim> vCornerCoords[],
751 const MathVector<refDim> vLocIP[],
752 const size_t nip,
753 LocalVector* u,
754 const MathMatrix<refDim, dim>* vJT = NULL) const
755 {
756 UG_ASSERT(dynamic_cast<elem_type*>(elem) != NULL, "Unsupported element type");
757 elem_type* element = static_cast<elem_type*>(elem);
758
759 // reference object id
761
762 const size_t numVertices = element->num_vertices();
763 // get domain of grid function
764 const domain_type& domain = *m_u->domain().get();
765
766 // get position accessor
767 typedef typename domain_type::position_accessor_type position_accessor_type;
768 const position_accessor_type& posAcc = domain.position_accessor();
769
770 position_accessor_type aaPos = m_u->domain()->position_accessor();
771
772 // coord and vertex array
775 DimCRFVGeometry<dim> crfvgeo;
776
778
779 for(size_t i = 0; i < numVertices; ++i){
780 vVrt[i] = element->vertex(i);
781 coCoord[i] = posAcc[vVrt[i]];
782 };
783
784 for (size_t i=0;i<domain_traits<dim>::MaxNumVerticesOfElem;i++)
785 grad[i]*=0.0;
786
787 // evaluate finite volume geometry
788 crfvgeo.update(elem, &(coCoord[0]), domain.subset_handler().get());
789
790 // Lagrange 1 trial space
791 const LocalShapeFunctionSet<dim>& lagrange1 =
792 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
793
794 std::vector<number> shapes;
795
796 for (size_t ip=0;ip<crfvgeo.num_scvf();ip++){
797 const typename DimCRFVGeometry<dim>::SCVF& scvf = crfvgeo.scvf(ip);
798 number pinter=0;
799 lagrange1.shapes(shapes,scvf.local_ip());
800 for (size_t sh=0;sh<numVertices;sh++)
801 pinter += m_p[vVrt[sh]]*shapes[sh];
802 for (int d=0;d<dim;d++){
803 number flux = pinter*scvf.normal()[d];
804 grad[scvf.from()][d]-=flux;
805 grad[scvf.to()][d]+=flux;
806 }
807 }
808
809 for (size_t i=0;i<crfvgeo.num_scv();i++){
810 grad[i]/=crfvgeo.scv(i).volume();
811 }
812
813 // evaluate source data
814 (*m_imSource)(vValue,
815 vGlobIP,
816 time, si,
817 elem,
818 vCornerCoords,
819 vLocIP,
820 nip,
821 u,
822 vJT);
823 for (size_t i=0;i<nip;i++)
824 vValue[i]+=grad[i];
825 }; // evaluate
826
827 void update(){
828 // get domain
829 domain_type& domain = *m_u->domain().get();
830 // create Multiindex
831 std::vector<DoFIndex> multInd;
833
834 // coord and vertex array
837
838 // get position accessor
839 typedef typename domain_type::position_accessor_type position_accessor_type;
840 const position_accessor_type& posAcc = domain.position_accessor();
841
842 // set volume and p values to zero
843 SetAttachmentValues(m_vol, m_u->template begin<Vertex>(), m_u->template end<Vertex>(), 0);
844 SetAttachmentValues(m_p, m_u->template begin<Vertex>(), m_u->template end<Vertex>(), 0);
845
846 // set p^{old} = p^{old} + p^{sep}
847 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si){
848 ElemIterator iter = m_u->template begin<elem_type>(si);
849 ElemIterator iterEnd = m_u->template end<elem_type>(si);
850 for( ;iter !=iterEnd; ++iter)
851 {
852 elem_type* elem = *iter;
853 m_u->inner_dof_indices(elem, _P_, multInd);
854 m_pOld[elem]+=DoFRef(*m_u,multInd[0]);
855 }
856 }
857 // compute pressure in vertices by averaging
858 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si){
859 ElemIterator iter = m_u->template begin<elem_type>(si);
860 ElemIterator iterEnd = m_u->template end<elem_type>(si);
861 for( ;iter !=iterEnd; ++iter)
862 {
863 elem_type* elem = *iter;
864 number pValue = m_pOld[elem];
865 const size_t numVertices = elem->num_vertices();
866 for(size_t i = 0; i < numVertices; ++i){
867 vVrt[i] = elem->vertex(i);
868 coCoord[i] = posAcc[vVrt[i]];
869 };
870 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
871 for(size_t i = 0; i < numVertices; ++i){
872 number scvVol = geo.scv(i).volume();
873 m_vol[vVrt[i]]+=scvVol;
874 m_p[vVrt[i]]+=scvVol*pValue;
875 }
876 }
877 }
878 PeriodicBoundaryManager* pbm = (domain.grid())->periodic_boundary_manager();
879 // go over all vertices and average
880 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si){
881 VertexIterator iter = m_u->template begin<Vertex>(si);
882 VertexIterator iterEnd = m_u->template end<Vertex>(si);
883 for( ;iter !=iterEnd; ++iter)
884 {
885 Vertex* vrt = *iter;
886 if (pbm && pbm->is_slave(vrt)) continue;
887 m_p[vrt]/=m_vol[vrt];
888 }
889 }
890 }
891
892 private:
893 static const size_t max_number_of_ips = 20;
894
895 public:
896 virtual void operator() (MathVector<dim>& value,
897 const MathVector<dim>& globIP,
898 number time, int si) const
899 {
900 UG_THROW("LevelSetUserData: Need element.");
901 }
902
903 virtual void operator() (MathVector<dim> vValue[],
904 const MathVector<dim> vGlobIP[],
905 number time, int si, const size_t nip) const
906 {
907 UG_THROW("LevelSetUserData: Need element.");
908 }
909
910 virtual void compute(LocalVector* u, GridObject* elem,
911 const MathVector<dim> vCornerCoords[], bool bDeriv = false)
912 {
913 const number t = this->time();
914 const int si = this->subset();
915 for(size_t s = 0; s < this->num_series(); ++s)
916 evaluate<dim>(this->values(s), this->ips(s), t, si,
917 elem, NULL, this->template local_ips<dim>(s),
918 this->num_ip(s), u);
919 }
920
922 virtual bool continuous() const {return false;}
923
925 virtual bool requires_grid_fct() const {return true;}
926};
927
928} // namespace NavierStokes
929} // end namespace ug
930
931
932#endif /* __H__UG__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__FVCR__PRESSURE_SEPARATION__ */
parameterString s
Definition Biogas.lua:2
function NavierStokes(fcts, subsets, discType)
T * get()
TData * values(size_t s)
TData & value(size_t s, size_t ip)
size_t num_ip(size_t s) const
size_t num_series() const
const MathVector< worldDim > & normal() const
const MathVector< dim > & local_ip() const
size_t num_scvf() const
size_t num_scv() const
const SCV & scv(size_t i) const
void update(GridObject *elem, const MathVector< worldDim > *vCornerCoords, const ISubsetHandler *ish=NULL)
const SCVF & scvf(size_t i) const
void update(GridObject *elem, const MathVector< worldDim > *vCornerCoords, const ISubsetHandler *ish=NULL)
const SCV & scv(size_t i) const
bool access(Grid &grid, TAttachment &a)
virtual ReferenceObjectID reference_object_id() const=0
int subset() const
size_t time_point()
number time() const
const MathVector< dim > * ips(size_t s) const
const MathVector< dim > & ip(size_t s, size_t ip) const
virtual void shapes(std::vector< std::vector< shape_type > > &vvShape, const std::vector< MathVector< dim > > &vLocPos) const=0
LocalVector & solution(size_t i)
std::size_t size() const
Definition pressure_separation.h:269
TGridFunction::template dim_traits< dim >::const_iterator ElemIterator
element iterator
Definition pressure_separation.h:297
domain_type::position_accessor_type position_accessor_type
position accessor type
Definition pressure_separation.h:277
aElementNumber m_pOld
Definition pressure_separation.h:305
virtual void compute(LocalVector *u, GridObject *elem, const MathVector< dim > vCornerCoords[], bool bDeriv=false)
Definition pressure_separation.h:574
virtual bool continuous() const
returns if provided data is continuous over geometric object boundaries
Definition pressure_separation.h:595
ANumber m_aVol
Definition pressure_separation.h:312
SeparatedPressureSource(SmartPtr< ApproximationSpace< domain_type > > approxSpace, SmartPtr< TGridFunction > spGridFct)
constructor
Definition pressure_separation.h:387
domain_type::grid_type grid_type
grid type
Definition pressure_separation.h:283
virtual ~SeparatedPressureSource()
Definition pressure_separation.h:407
PeriodicAttachmentAccessor< Vertex, ANumber > aVertexNumber
MathVector<dim> attachment.
Definition pressure_separation.h:293
virtual void compute(LocalVectorTimeSeries *u, GridObject *elem, const MathVector< dim > vCornerCoords[], bool bDeriv=false)
Definition pressure_separation.h:584
ANumber m_aP
Definition pressure_separation.h:308
TGridFunction::algebra_type algebra_type
algebra type
Definition pressure_separation.h:274
TGridFunction::template traits< Vertex >::const_iterator VertexIterator
vertex iterator
Definition pressure_separation.h:300
static const size_t max_number_of_ips
Definition pressure_separation.h:557
SmartPtr< TGridFunction > m_u
Definition pressure_separation.h:316
void set_source(SmartPtr< CplUserData< MathVector< dim >, dim > > data)
Definition pressure_separation.h:335
SmartPtr< ApproximationSpace< domain_type > > m_spApproxSpace
Definition pressure_separation.h:319
Grid::AttachmentAccessor< elem_type, ANumber > aElementNumber
Definition pressure_separation.h:294
aVertexNumber m_vol
Definition pressure_separation.h:313
TGridFunction::template dim_traits< dim >::grid_base_object elem_type
element type
Definition pressure_separation.h:286
aVertexNumber m_p
Definition pressure_separation.h:309
static const int dim
world dimension
Definition pressure_separation.h:280
void set_source(number f_x, number f_y, number f_z)
Definition pressure_separation.h:362
virtual bool requires_grid_fct() const
returns if grid function is needed for evaluation
Definition pressure_separation.h:598
TGridFunction::domain_type domain_type
domain type
Definition pressure_separation.h:271
SmartPtr< CplUserData< MathVector< dim >, dim > > m_imSource
Data import for source.
Definition pressure_separation.h:330
void evaluate(MathVector< dim > vValue[], const MathVector< dim > vGlobIP[], number time, int si, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< refDim > vLocIP[], const size_t nip, LocalVector *u, const MathMatrix< refDim, dim > *vJT=NULL) const
Definition pressure_separation.h:410
grid_type * m_grid
Definition pressure_separation.h:322
void update()
Definition pressure_separation.h:491
virtual void operator()(MathVector< dim > &value, const MathVector< dim > &globIP, number time, int si) const
Definition pressure_separation.h:560
void set_source(number f_x)
Definition pressure_separation.h:340
ANumber m_aPOld
Definition pressure_separation.h:304
void set_source(number f_x, number f_y)
Definition pressure_separation.h:349
static const size_t _P_
Definition pressure_separation.h:325
Definition pressure_separation.h:605
virtual bool continuous() const
returns if provided data is continuous over geometric object boundaries
Definition pressure_separation.h:922
SmartPtr< TGridFunction > m_u
Definition pressure_separation.h:652
Grid::AttachmentAccessor< elem_type, ANumber > aElementNumber
Definition pressure_separation.h:630
ANumber m_aP
Definition pressure_separation.h:644
virtual bool requires_grid_fct() const
returns if grid function is needed for evaluation
Definition pressure_separation.h:925
aVertexNumber m_vol
Definition pressure_separation.h:649
void set_source(number f_x, number f_y)
Definition pressure_separation.h:685
SmartPtr< CplUserData< MathVector< dim >, dim > > m_imSource
Data import for source.
Definition pressure_separation.h:666
TGridFunction::template traits< Vertex >::const_iterator VertexIterator
vertex iterator
Definition pressure_separation.h:636
domain_type::grid_type grid_type
grid type
Definition pressure_separation.h:619
void update()
Definition pressure_separation.h:827
TGridFunction::domain_type domain_type
domain type
Definition pressure_separation.h:607
void set_source(number f_x)
Definition pressure_separation.h:676
virtual void compute(LocalVector *u, GridObject *elem, const MathVector< dim > vCornerCoords[], bool bDeriv=false)
Definition pressure_separation.h:910
virtual ~SeparatedPressureSourceInter()
Definition pressure_separation.h:743
aVertexNumber m_p
Definition pressure_separation.h:645
void evaluate(MathVector< dim > vValue[], const MathVector< dim > vGlobIP[], number time, int si, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< refDim > vLocIP[], const size_t nip, LocalVector *u, const MathMatrix< refDim, dim > *vJT=NULL) const
Definition pressure_separation.h:746
TGridFunction::algebra_type algebra_type
algebra type
Definition pressure_separation.h:610
aElementNumber m_pOld
Definition pressure_separation.h:641
ANumber m_aPOld
Definition pressure_separation.h:640
TGridFunction::template dim_traits< dim >::const_iterator ElemIterator
element iterator
Definition pressure_separation.h:633
void set_source(SmartPtr< CplUserData< MathVector< dim >, dim > > data)
Definition pressure_separation.h:671
domain_type::position_accessor_type position_accessor_type
position accessor type
Definition pressure_separation.h:613
grid_type * m_grid
Definition pressure_separation.h:658
SmartPtr< ApproximationSpace< domain_type > > m_spApproxSpace
Definition pressure_separation.h:655
TGridFunction::template dim_traits< dim >::grid_base_object elem_type
element type
Definition pressure_separation.h:622
SeparatedPressureSourceInter(SmartPtr< ApproximationSpace< domain_type > > approxSpace, SmartPtr< TGridFunction > spGridFct)
constructor
Definition pressure_separation.h:723
PeriodicAttachmentAccessor< Vertex, ANumber > aVertexNumber
MathVector<dim> attachment.
Definition pressure_separation.h:629
ANumber m_aVol
Definition pressure_separation.h:648
void set_source(number f_x, number f_y, number f_z)
Definition pressure_separation.h:698
bool access(Grid &g, TAttachment &a)
bool is_slave(TElem *) const
SmartPtr< TGrid > grid()
size_t num_subsets() const
void SetAttachmentValues(TAttachmentAccessor &aaVal, TIter elemsBegin, TIter elemsEnd, const TVal &val)
#define UG_ASSERT(expr, msg)
#define UG_THROW(msg)
double number
bool solveLS(std::vector< number > &x, const std::vector< number > &matField, const std::vector< number > &b)
Definition pressure_separation.h:82
int bubblesort(std::vector< int > &list, std::vector< number > array)
Definition pressure_separation.h:238
bool multMatVec(const std::vector< number > &avec, const std::vector< number > &b, std::vector< number > &c, size_t m, size_t n)
Definition pressure_separation.h:67
void computeElemBarycenter(MathVector< dim > &bary, TElem *elem, TPos posA)
Definition pressure_separation.h:58
bool leastSquares(std::vector< number > &x, const std::vector< number > &mField, const std::vector< number > &b)
Definition pressure_separation.h:213
ReferenceObjectID
const number & DoFRef(const TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
bool resize(size_t newRows, size_t newCols)