Plugins
Loading...
Searching...
No Matches
level_set_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2011-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
34#ifndef LEVEL_SET_UTIL_IMPL_H_
35#define LEVEL_SET_UTIL_IMPL_H_
36
37#include "lib_disc/local_finite_element/local_finite_element_provider.h"
38#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
39#include "lib_disc/reference_element/reference_element.h"
40#include "lib_grid/algorithms/attachment_util.h"
41
42#include <algorithm>
43
44namespace ug{
45namespace LevelSet{
46
51template<typename TGridFunction>
52bool FV1LevelSetDisc<TGridFunction>::limit_grad
53(
54 TGridFunction& uOld,
55 aaGrad& aaGrad
56)
57{
58 // get grid
59 grid_type& grid = *uOld.domain()->grid();
60
61 position_accessor_type& aaPos = uOld.domain()->position_accessor();
62
63 std::vector<DoFIndex> ind;
64
65 // create Attachment for scv-volume size
66 ANumber aMax;
67 ANumber aMin;
68
69 // attach to grid
70 grid.attach_to_vertices(aMin);
71 grid.attach_to_vertices(aMax);
72
73 // get attachment accessor to access values
74 Grid::VertexAttachmentAccessor<ANumber> aaMin(grid, aMin);
75 Grid::VertexAttachmentAccessor<ANumber> aaMax(grid, aMax);
76 for (int si=0;si<uOld.num_subsets();++si)
77 {
78 for(VertexConstIterator iter = uOld.template begin<Vertex>(si);
79 iter !=uOld.template end<Vertex>(si); ++iter)
80 {
81 Vertex* vrt = *iter;
82 MathVector<dim> coord;
83 coord = aaPos[vrt];
84 // read indices on vertex
85 // get vector holding all indices on the vertex
86 uOld.inner_dof_indices(vrt, 0, ind);
87 aaMax[vrt] = DoFRef(uOld, ind[0]);
88 aaMin[vrt] = DoFRef(uOld, ind[0]);
89 }
90 }
91 for (int si=0;si<uOld.num_subsets();++si)
92 {
93 //UG_LOG("si " << si << "\n");
94 for(EdgeConstIterator iter = uOld.template begin<Edge>(si) ;
95 iter !=uOld.template end<Edge>(si); ++iter)
96 {
97 Edge* edge = *iter;
98 Vertex* vi=edge->vertex(0);
99 Vertex* vj=edge->vertex(1);
100 uOld.inner_dof_indices(vi, 0, ind);
101 number ui = DoFRef(uOld, ind[0]);
102 uOld.inner_dof_indices(vj, 0, ind);
103 number uj = DoFRef(uOld, ind[0]);
104 //UG_LOG("edge " << aaPos[vi] << "-" << aaPos[vj] << " [" << ui << " " << uj << "]\n");
105 if (uj<aaMin[vi])
106 aaMin[vi]=uj;
107 if (uj>aaMax[vi])
108 aaMax[vi]=uj;
109 if (ui<aaMin[vj])
110 aaMin[vj]=ui;
111 if (ui>aaMax[vj])
112 aaMax[vj]=ui;
113 }
114
115 };
116/*
117 for (int si=0;si<2;++si)
118 {
119 for(VertexConstIterator iter = uOld.template begin<Vertex>(si) ;
120 iter !=uOld.template end<Vertex>(si); ++iter)
121 {
122 Vertex* vrt = *iter;
123 MathVector<dim> coord;
124 coord = aaPos[vrt];
125 // read indices on vertex
126 // get vector holding all indices on the vertex
127 uOld.inner_dof_indices(vrt, 0, ind);
128 //UG_LOG(coord << " min=" << aaMin[vrt] << " max=" << aaMax[vrt] << "\n");
129 }
130 }
131*/
132 for (int si=0;si<1;++si)
133 {
134 for(EdgeConstIterator iter = uOld.template begin<Edge>(si) ;
135 iter !=uOld.template end<Edge>(si); ++iter)
136 {
137 Edge* edge = *iter;
138 Vertex* vi=edge->vertex(0);
139 Vertex* vj=edge->vertex(1);
140 MathVector<dim> coordi,coordj,coordij,distVec,gradi,gradj;
141 gradi = aaGrad[vi];
142 gradj = aaGrad[vj];
143 coordi = aaPos[vi];
144 coordj = aaPos[vj];
145 uOld.inner_dof_indices(vi, 0, ind);
146 number ui = DoFRef(uOld, ind[0]);
147 uOld.inner_dof_indices(vj, 0, ind);
148 number uj = DoFRef(uOld, ind[0]);
149 VecScaleAdd(coordij,0.5,coordi,0.5,coordj);
150 VecSubtract(distVec, coordij,coordi);
151 number uij = ui + distVec*gradi;
152 number alpha = 1;
153 if (uij>ui)
154 {
155 if (uij>aaMax[vi]) alpha=(aaMax[vi]-ui)/(distVec*gradi);
156 if (alpha<1)
157 {
158 //UG_LOG("edge " << coordi << " " << coordj << "\n");
159 //UG_LOG(coordi << " u " << ui << " uij " << uij << " max " << aaMax[vi] << " alpha " << alpha << "\n");
160 aaGrad[vi]*=alpha;
161 };
162 }
163 else
164 {
165 if (uij<aaMin[vi]) alpha=(aaMin[vi]-ui)/(distVec*gradi);
166 if (alpha<1)
167 {
168 //UG_LOG("edge " << coordi << " " << coordj << "\n");
169 //UG_LOG(coordi << " u " << ui << " uij " << uij << " min " << aaMax[vi] << " alpha " << alpha << "\n");
170 aaGrad[vi]*=alpha;
171 };
172 };
173 VecSubtract(distVec, coordij,coordj);
174 uij = uj + distVec*gradj;
175 alpha = 1;
176 if (uij>uj)
177 {
178 if (uij>aaMax[vj]) alpha=(aaMax[vj]-uj)/(distVec*gradj);
179 if (alpha<1)
180 {
181 //UG_LOG("-- edge " << coordi << " " << coordj << "\n");
182 //UG_LOG(coordj << " u " << uj << " uij " << uij << " max " << aaMax[vj] << " alpha " << alpha << "\n");
183 aaGrad[vj]*=alpha;
184 };
185 }
186 else
187 {
188 if (uij<aaMin[vj]) alpha=(aaMin[vj]-uj)/(distVec*gradj);
189 if (alpha<1)
190 {
191 //UG_LOG("-- edge " << coordi << " " << coordj << "\n");
192 //UG_LOG(coordj << " u " << uj << " uij " << uij << " min " << aaMax[vj] << " alpha " << alpha << "\n");
193 aaGrad[vj]*=alpha;
194 };
195 };
196 // UG_LOG(" coord vertex 0 " << aaPos[v0] << " coord vertex 1 " << aaPos[v1] << "\n");
197 }
198
199 Vertex* vVrt[domain_traits<dim>::MaxNumVerticesOfElem];
200 // coord and vertex array
201 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
202 // coord and vertex array
203 MathVector<dim> grad[domain_traits<dim>::MaxNumVerticesOfElem];
204 // values of the grid function
205 number u[domain_traits<dim>::MaxNumVerticesOfElem];
206 // get iterators
207 ElemIterator iter = uOld.template begin<ElemType>(si);
208 ElemIterator iterEnd = uOld.template end<ElemType>(si);
209 // loop elements of dimension
210 for( ;iter !=iterEnd; ++iter)
211 {
212 // get Elem
213 ElemType* elem = *iter;
214
215 // get position accessor
216 const position_accessor_type& aaPos = uOld.domain()->position_accessor();
217
218 // compute center of mass
219 MathVector<dim> center;
220 center=0;
221 size_t noc=elem->num_vertices();
222 for(size_t i = 0; i < noc; ++i)
223 {
224 vVrt[i] = elem->vertex(i);
225 coCoord[i] = aaPos[vVrt[i]];
226 grad[i] = aaGrad[vVrt[i]] ;
227 VecAppend(center,coCoord[i]);
228 uOld.inner_dof_indices(vVrt[i], 0, ind);
229 u[i]=DoFRef(uOld, ind[0]);
230 };
231 center/=noc;
232 for (size_t i=0;i<noc;++i)
233 {
234 number alpha=1;
235 MathVector<dim> distVec;
236 number uCenter;
237 VecSubtract(distVec,center,coCoord[i]);
238 uCenter = u[i] + distVec*grad[i];
239 if (uCenter>u[i])
240 {
241 if (uCenter>aaMax[vVrt[i]]) alpha=(aaMax[vVrt[i]]-u[i])/(distVec*grad[i]);
242 if (alpha<1)
243 {
244 // UG_LOG("* " << coCoord[i] << " uCenter " << uCenter << " ui " << u[i] << " max " << aaMax[vVrt[i]] << " alpha " << alpha << "\n");
245 aaGrad[vVrt[i]]*=alpha;
246 };
247 }
248 else
249 {
250 if (uCenter<aaMin[vVrt[i]]) alpha=(aaMin[vVrt[i]]-u[i])/(distVec*grad[i]);
251 if (alpha<1)
252 {
253 // UG_LOG("*#* " << coCoord[i] << " uCenter " << uCenter << " ui " << u[i] << " min " << aaMin[vVrt[i]] << " alpha " << alpha << "\n");
254 aaGrad[vVrt[i]]*=alpha;
255 };
256 };
257 };
258 };
259 };
260 // detach from grid
261 grid.detach_from_vertices(aMin);
262 grid.detach_from_vertices(aMax);
263 return true;
264};
265
266
273template<typename TGridFunction>
274template <typename TElem>
275bool FV1LevelSetDisc<TGridFunction>::assemble_element
276(
277 TElem& elem,
278 DimFV1Geometry<dim>& geo,
279 grid_type& grid,
280 TGridFunction& uNew,
281 const TGridFunction& uOld,
282 aaGrad& aaGradient,
283 aaVol& aaVolume
284)
285{
286// a large enough number for max. number of corners and subcontrol volume faces
287 static const size_t maxNumCo = 20;
288
289// get domain
290 domain_type& domain = *uNew.domain().get();
291
292// create Multiindex
293 std::vector<DoFIndex> multInd;
294
295 // hard code function (fct=0)
296 //\todo: generalize
297// size_t fct=0;
298
299// get position accessor
300 const position_accessor_type& aaPos = domain.position_accessor();
301
302// coord and vertex array
303 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
304 Vertex* vVrt[domain_traits<dim>::MaxNumVerticesOfElem];
305
306// get vertices and extract corner coordinates
307 const size_t numVertices = elem->num_vertices();
308 for(size_t i = 0; i < numVertices; ++i)
309 {
310 vVrt[i] = elem->vertex(i);
311 coCoord[i] = aaPos[vVrt[i]];
312 };
313
314// update fv geometry
315 geo.update(elem, &(coCoord[0]), uOld.domain()->subset_handler().get());
316
317// fill node value vector
318 std::vector<number> uValue(geo.num_scv());
319 size_t noc = geo.num_scv();
320 for (size_t i=0; i < noc; i++)
321 {
322 // if (dd.template inner_dof_indices<Vertex>(vVrt[i], 0, multInd) != 1) return false;
323 uOld.inner_dof_indices(vVrt[i], 0, multInd);
324 uValue[i]=DoFRef(uOld, multInd[0]);
325 }
326// fill grad vector
327 MathVector<dim> grad[maxNumCo];
328 for (size_t i=0; i < noc; i++)
329 {
330 // if (dd.template inner_dof_indices<Vertex>(vVrt[i], 0, multInd) != 1) return false;
331 uOld.inner_dof_indices(vVrt[i], 0, multInd);
332 grad[i]= aaGradient[vVrt[i]];
333 }
334
335// fill corner velocity vector
336 MathVector<dim> coVelocity[maxNumCo];
337 for (size_t i=0; i<noc; ++i)
338 coVelocity[i]=0;
339 if (m_gamma!=0)
340 {
341 const int si = 0; //TODO this should be corrected
342 if (m_imVelocity->requires_grid_fct())
343 {
344 // create storage
345 LocalIndices localind;
346 LocalVector localu;
347
348 // get global indices
349 uNew.indices(elem, localind);
350
351 // adapt local algebra
352 localu.resize(localind);
353
354 // read local values of u
355 GetLocalVector(localu, uNew);
356
357 (*m_imVelocity)(coVelocity, geo.scv_global_ips(), m_time, si,
358 elem,coCoord, geo.scv_local_ips(),geo.num_sh(), &localu);
359 }
360 else
361 {
362 // see user_data.h : 410
363 (*m_imVelocity)(coVelocity, geo.scv_global_ips(), m_time, si, geo.num_sh());
364 }
365 if (m_gamma!=1) for (size_t i=0; i < noc; i++) coVelocity[i]*=m_gamma;
366 }
367 if (m_delta!=0)
368 {
369 for (size_t i=0; i < noc; i++)
370 {
371 number vnorm = VecLength(grad[i]);
372 if (vnorm>1e-15) for (int j=0; j < dim; j++) coVelocity[i][j] += m_delta/vnorm*grad[i][j];
373 };
374 };
375
376// compute the ip velocity from the corner velocity
377 MathVector<dim> ipVelocity[maxNumCo];
378 for (size_t ip=0;ip < geo.num_scvf();ip++)
379 {
380 ipVelocity[ip] = 0;
381 const typename DimFV1Geometry<dim>::SCVF& scvf = geo.scvf(ip);
382 for (size_t co=0; co < noc; co++)
383 {
384 for (int j=0; j<dim; j++)
385 ipVelocity[ip][j] += scvf.shape(co)*coVelocity[co][j];
386 };
387 }
388
389// fill source vector
390 std::vector<number> coSource(noc);
391 const int si = 0; //TODO this should be corrected
392 (*m_imSource)(&coSource[0], geo.scv_global_ips(), m_time, si, geo.num_sh());
393
394// compute fluxes
395 size_t base;
396 number flux;
397 MathVector<dim> distVec;
398 for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
399 {
400 MathVector<dim> bNode;
401 // get current SCVF
402 const typename DimFV1Geometry<dim>::SCVF& scvf = geo.scvf(ip);
403 MathVector<dim> ipCoord = scvf.global_ip();
404 size_t from = scvf.from();
405 size_t to = scvf.to();
406 if (scvf.normal()*ipVelocity[ip]>0)
407 base = from;
408 else
409 base = to;
410 VecSubtract(distVec, ipCoord,coCoord[base]);
411 // flux = v * n * u_{ip(i)}^{n+0.5}
412 flux = m_dt*(ipVelocity[ip]*scvf.normal())*( uValue[base] + (distVec*grad[base]) + 0.5*m_dt*(coSource[base] - (grad[base]*coVelocity[base])) );
413 //UG_LOG(ip << " flux=" << flux << "\n");
414 uOld.inner_dof_indices(vVrt[from], 0, multInd);
415 DoFRef(uNew, multInd[0])-=flux/aaVolume[ vVrt[from] ];
416 if (! m_divFree)
417 {
418 DoFRef(uNew, multInd[0])+=m_dt*(ipVelocity[ip]*scvf.normal())*(uValue[from] + 0.5*m_dt*(coSource[from] - (grad[from]*coVelocity[from])))/aaVolume[ vVrt[from] ];
419 //DofRef(uNew[multInd[0][0]],multInd[0][1])+=m_dt*(ipVelocity[ip]*scvf.normal())*(uValue[from] + 0.0*m_dt*(coSource[from] - (grad[from]*coVelocity[from])))/aaVolume[ vVrt[from] ];
420 };
421 uOld.inner_dof_indices(vVrt[to], 0, multInd);
422 DoFRef(uNew, multInd[0])+=flux/aaVolume[ vVrt[to] ];
423 if (! m_divFree)
424 {
425 DoFRef(uNew, multInd[0])-=m_dt*(ipVelocity[ip]*scvf.normal())*(uValue[to] + 0.5*m_dt*(coSource[to] - (grad[to]*coVelocity[to])))/aaVolume[ vVrt[to] ];
426 };
427 number localCFL = std::max(m_dt*std::abs(ipVelocity[ip]*scvf.normal())/aaVolume[ vVrt[from] ],m_dt*std::abs(ipVelocity[ip]*scvf.normal())/aaVolume[ vVrt[to] ] );
428 if (localCFL>m_maxCFL)
429 {
430 m_maxCFL = localCFL;
431 };
432 };
433
434// boundary
435 if (geo.num_bf()>0)
436 {
437 for (int si=0;si<uNew.num_subsets();++si)
438 {
439 //UG_LOG("si=" << si << " num_bf(si)=" << geo.num_bf(si) << "\n");
440 for(size_t i = 0; i < geo.num_bf(si); ++i)
441 {
442 // get current BF
443 const typename DimFV1Geometry<dim>::BF& bf = geo.bf(si, i);
444 const size_t nodeID = bf.node_id();
445 MathVector<dim> bipVelocity;
446 number bipU=0;
447 number bipSource=0;
448 MathVector<dim> bipGrad;
449 bipVelocity=0;
450 bipGrad=0;
451 const MathVector<dim>* globalGradVec = bf.global_grad_vector();
452 bipVelocity=0;
453 for (size_t co=0;co<noc;co++)
454 {
455 //UG_LOG("corner " << co << " num_sh " << bf.num_sh() << "\n");
456 for (int j=0;j<dim;j++)
457 {
458 bipVelocity[j] += bf.shape(co) * coVelocity[co][j];
459 bipGrad[j] += bf.shape(co) * globalGradVec[co][j];
460 };
461 bipU += bf.shape(co) * uValue[co];
462 bipSource += bf.shape(co) * coSource[co];
463 }
464 //flux = m_dt*(bipVelocity*bf.normal())*( bipU + 0.5*m_dt*(bipSource - (bipGrad*bipVelocity)) );
465 flux = m_dt*(bipVelocity*bf.normal())*uValue[nodeID];// first order approximation
466 uOld.inner_dof_indices(vVrt[nodeID], 0, multInd);
467 DoFRef(uNew, multInd[0])-=flux/aaVolume[ vVrt[nodeID] ];
468 if (!m_divFree)
469 {
470 DoFRef(uNew, multInd[0])+=m_dt*(bipVelocity*bf.normal())*(uValue[nodeID] )/aaVolume[ vVrt[nodeID] ];// first order approximation
471 };
472 };
473 };
474 };
475
476 return true;
477}
478
483template<typename TGridFunction>
484bool FV1LevelSetDisc<TGridFunction>::calculate_vertex_grad_vol
485(
486 TGridFunction& u,
487 aaGrad& aaGradient,
488 aaVol& aaVolume
489)
490{
491// get domain
492 domain_type& domain = *u.domain().get();
493
494// get grid of domain
495 typename domain_type::grid_type& grid = *domain.grid();
496
497// create Multiindex
498 std::vector<DoFIndex> multInd;
499
500// create a FV Geometry for the dimension
501 DimFV1Geometry<dim> geo;
502
503// hard code function (fct=0)
504//\todo: generalize
505 size_t fct=0;
506
507 // initialize attachment value
508 SetAttachmentValues(aaVolume, grid.vertices_begin(), grid.vertices_end(), 0);
509 SetAttachmentValues(aaGradient, grid.vertices_begin(), grid.vertices_end(), 0);
510
511 // coord and vertex array
512 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
513 Vertex* vVrt[domain_traits<dim>::MaxNumVerticesOfElem];
514
515// sum up all contributions of the sub control volumes to one vertex in an attachment
516 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
517 {
518 // skip boundary
519 if (m_dirichlet_sg.size()!=0) if (m_dirichlet_sg.contains(si)) continue;
520 if (m_neumann_sg.size()!=0) if (m_neumann_sg.contains(si)) continue;
521 if (m_inactive_sg.size()!=0) if (m_inactive_sg.contains(si)) continue;
522
523 // get iterators
524 ElemIterator iter = u.template begin<ElemType>(si);
525 ElemIterator iterEnd = u.template end<ElemType>(si);
526
527 // loop elements of dimension
528 for( ;iter !=iterEnd; ++iter)
529 {
530 // get Elem
531 ElemType* elem = *iter;
532
533 // get position accessor
534 const position_accessor_type& aaPos = domain.position_accessor();
535
536 // get vertices and extract corner coordinates
537 const size_t numVertices = elem->num_vertices();
538 for(size_t i = 0; i < numVertices; ++i)
539 {
540 vVrt[i] = elem->vertex(i);
541 coCoord[i] = aaPos[vVrt[i]];
542 };
543
544 // evaluate finite volume geometry
545 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
546
547 //UG_LOG("Num Verts loaded: "<<vVrt.size()<<"\n");
548 //UG_LOG("Num SCV computed: "<<geo.num_scv()<<"\n");
549
550 number uValue[domain_traits<dim>::MaxNumVerticesOfElem];
551
552 // read indices on vertex
553 size_t noc = geo.num_scv();
554 for (size_t i=0;i < noc;i++)
555 {
556 // get indices of function fct on vertex
557 u.inner_dof_indices(vVrt[i], fct, multInd);
558
559 // read value of index from vector
560 uValue[i]=DoFRef(u, multInd[0]);
561
562 // debug log
563 //UG_LOG("corner " << i << " " << uValue[i] << "\n");
564 }
565
566 // storage for global gradient
567 MathVector<dim> globalGrad;
568
569 // loop corners
570 for (size_t i=0;i < noc;i++)
571 {
572 // get scv for sh
573 const typename DimFV1Geometry<dim>::SCV& scv = geo.scv(i);
574
575 // debug log
576 //UG_LOG("gradient for corner " << i << "\n");
577
578 // reset global gradient
579 globalGrad = 0.0;
580
581 // sum up gradients of shape functions in corner
582 for(size_t sh = 0 ; sh < noc; ++sh)
583 {
584 //UG_LOG("local grad " << sh << " : " << scv.local_grad(sh) << "\n");
585 //UG_LOG("unscaled global grad " << sh << " = " << scv.global_grad(sh) << "\n");
586 //UG_LOG("uvalue(" << sh << ") =" << uValue[sh] << "\n");
587
588 VecScaleAppend(globalGrad, uValue[sh], scv.global_grad(sh));
589 }
590
591 // volume of scv
592 number vol = scv.volume();
593
594 //UG_LOG("*** global grad " << i << ": " << globalGrad << "\n");
595
596 // scale gradient by volume
597 globalGrad *= vol;
598
599 // add both values to attachements
600 aaGradient[vVrt[i]] += globalGrad;
601 aaVolume[vVrt[i]] += vol;
602 };
603 }
604 }
605
606// divide the gradients by the volumes
607 // int count=0;
608 // position_accessor_type aaPos = u.domain()->position_accessor();
609 for (int si=0;si < u.num_subsets();++si)
610 {
611 //UG_LOG("si " << si << "\n");
612 for(VertexConstIterator iter = u.template begin<Vertex>(si);
613 iter != u.template end<Vertex>(si); ++iter)
614 {
615 // get vertex
616 Vertex* vrt = *iter;
617 if (aaVolume[vrt]!=0)
618 {
619 (aaGradient[vrt]) /= aaVolume[vrt];
620 };
621 //exact[0] = cos(coord[0]);// 6*coord[0];
622 //exact[1] = -4*sin(coord[1]);//-4*coord[1];
623 //number gError = sqrt( (exact[0]-aaGradient[vrt][0])*(exact[0]-aaGradient[vrt][0]) + (exact[1]-aaGradient[vrt][1])*(exact[1]-aaGradient[vrt][1]) );
624 //UG_LOG(count << "[ " << coord[0] << "," << coord[1] << " ] vol= " << aaVolume[vrt] << " " << "grad= ["
625 // << aaGradient[vrt][0] << "," << aaGradient[vrt][1] << "] exact grad = ["
626 // << exact[0] << "," << exact[1] << "] error: " << gError << "\n");
627 // count++;
628 }
629 }
630
631 //UG_LOG("#*#*#*#\n");
632 return true ;
633}
634
638template<typename TGridFunction>
639bool FV1LevelSetDisc<TGridFunction>::assign_dirichlet
640(
641 TGridFunction& numsol
642)
643{
644// get domain of grid function
645 domain_type& domain = *numsol.domain().get();
646
647// UG_LOG("dirichlet\n");
648
649 // UG_LOG("nr dir ss " << m_dirichlet_sg.size() << "\n");
650
651 for(size_t i = 0; i < m_dirichlet_sg.size(); ++i)
652 {
653 const int si = m_dirichlet_sg[i];
654 // UG_LOG("Dirichlet boundary is: "<<si<< "\n");
655 for(VertexConstIterator iter = numsol.template begin<Vertex>(si);
656 iter != numsol.template end<Vertex>(si); ++iter)
657 {
658 // get vertex
659 Vertex* vrt = *iter;
660 number exactVal;
661 position_accessor_type aaPos = domain.position_accessor();
662
663 // get vector holding all indices on the vertex
664 std::vector<DoFIndex> ind;
665
666 const size_t numInd = numsol.inner_dof_indices(vrt, 0, ind);
667
668 // check indices
669 if(numInd != 1) {UG_LOG("ERROR: Wrong number of indices!"); return false;}
670
671 (*m_imDirichlet)(&exactVal,&aaPos[vrt],m_time,si,1);
672 DoFRef(numsol, ind[0]) = exactVal;
673
674 // MathVector<dim> coord = aaPos[vrt];
675 // UG_LOG("coord " << coord[0] << "," << coord[1] << " <> " << BlockRef(numsol[ind[0][0]],ind[0][1]) << "\n");
676
677 //if ((coord[0]==-1)||(coord[0]==1)||(coord[1]==-1)||(coord[1]==1)){
678 //BlockRef(numsol[ind[0][0]],ind[0][1]) = exactVal;
679 //};
680 //if ((coord[0]==0)||(coord[0]==1)||(coord[1]==0)||(coord[1]==1)){
681 //};
682 }
683 }
684
685 return true;
686}
687
694template<typename TGridFunction>
695bool FV1LevelSetDisc<TGridFunction>::advect_lsf
696(
697 TGridFunction& uNew,
698 TGridFunction& uOld
699)
700{
701// get domain of grid function
702 domain_type& domain = *uNew.domain().get();
703
704// get grid of domain
705 grid_type& grid = *domain.grid();
706
707// create Attachment for scv-volume size
708 ANumber aScvVolume;
709
710// typedef of gradient attachment
711 typedef Attachment<MathVector<dim> > AGradient;
712
713// create Attachment for gradient
714 AGradient aGradient;
715
716// attach to grid
717 grid.attach_to_vertices(aScvVolume);
718 grid.attach_to_vertices(aGradient);
719
720// get attachment accessor to access values
721 Grid::VertexAttachmentAccessor<ANumber> aaVolume(grid, aScvVolume);
722 Grid::VertexAttachmentAccessor<AGradient> aaGradient(grid, aGradient);
723
724// initialize attachment value
725 SetAttachmentValues(aaVolume, grid.vertices_begin(), grid.vertices_end(), 0);
726 SetAttachmentValues(aaGradient, grid.vertices_begin(), grid.vertices_end(), 0);
727
728
729// the CFL constant to compute
730 m_maxCFL=0;
731
732 //UG_LOG("***************************************************\n");
733 //UG_LOG("***************************************************\n");
734 //UG_LOG("***************************************************\n");
735
736 VecAssign(uNew,uOld);
737 MathVector<dim> coord;
738 std::vector<DoFIndex> ind;
739 position_accessor_type aaPos = domain.position_accessor();
740
741// compute time steps
742 for (size_t step=0; step<m_nrOfSteps; step++)
743 {
744 // calculate scv volume and the gradient
745 if (! calculate_vertex_grad_vol(uNew,aaGradient, aaVolume)) {UG_LOG("ERROR: gradient computation failed!"); };
746 // SetAttachmentValues(aaGradient, grid.vertices_begin(), grid.vertices_end(), 0); // for debug set gradient to 0
747 if (m_limiter)
748 limit_grad(uNew,aaGradient);
749
750 // UG_LOG("num_subsets: " << uOld.num_subsets() << "\n");
751
752 // loop over subsets to compute the new solution
753 for (int si=0;si<uOld.num_subsets();++si)
754 {
755 // skip boundaries
756 //UG_LOG("si " << si << "\n");
757 if (m_dirichlet_sg.size()!=0) if (m_dirichlet_sg.contains(si)) continue;
758 if (m_neumann_sg.size()!=0) if (m_neumann_sg.contains(si)) continue;
759 if (m_inactive_sg.size()!=0) if (m_inactive_sg.contains(si)) continue;
760
761 //UG_LOG("... \n");
762 //UG_LOG("***************************************************\n");
763 //UG_LOG("***********************" << si << "**************************\n");
764 //UG_LOG("***************************************************\n");
765
766 // get iterators
767 ElemIterator iter = uNew.template begin<ElemType>(si);
768 ElemIterator iterEnd = uNew.template end<ElemType>(si);
769
770 DimFV1Geometry<dim> geo;
771
772 // flag given neumann bnd subsets at the geometry, such that the
773 // geometry produces boundaryfaces (BF) for all sides of the
774 // element, that is in one of the subsets
775 for(size_t i = 0; i < m_neumann_sg.size(); ++i)
776 {
777 const int bndSi = m_neumann_sg[i];
778// UG_LOG("Neumann boundary is: "<<bndSi<< "\n");
779 geo.add_boundary_subset(bndSi);
780 }
781
782 // loop elements compute the new solution
783 for( ;iter !=iterEnd; ++iter)
784 {
785 // get Elem
786 //UG_LOG("*** ELEM ***\n");
787 ElemType* elem = *iter;
788 //UG_LOG("element \n");
789 // uNew = uOld
790 assemble_element(elem, geo, grid, uNew, uOld, aaGradient, aaVolume);
791 };
792 };
793
794 // take into account the source at the vertices
795 for (int si=0;si<uOld.num_subsets();++si)
796 {
797 if (m_dirichlet_sg.size()!=0) if (m_dirichlet_sg.contains(si)) continue;
798 if (m_inactive_sg.size()!=0) if (m_inactive_sg.contains(si)) continue;
799 std::vector<number> sourceValue(1);
800 MathVector<dim> sourceCo[1];
801 for(VertexConstIterator iter = uNew.template begin<Vertex>(si);
802 iter != uNew.template end<Vertex>(si); ++iter)
803 {
804 Vertex* vrt = *iter;
805 sourceCo[0]= aaPos[vrt];
806 (*m_imSource)(&sourceValue[0],sourceCo,m_time,si,1);
807 uNew.inner_dof_indices(vrt, 0, ind);
808 DoFRef(uNew, ind[0]) += m_dt*sourceValue[0];
809 };
810 }
811
812 // the new solution computed
813 m_time += m_dt;
814 m_timestep_nr++;
815
816 // set the Dirichlet values
817 assign_dirichlet(uNew);
818
819 // overwrite inactive nodes with old solution
820 for(size_t i = 0; i < m_inactive_sg.size(); ++i)
821 {
822 const int si = m_inactive_sg[i];
823 UG_LOG("inactive si: " << si << "\n");
824 for(VertexConstIterator iter = uNew.template begin<Vertex>(si);
825 iter != uNew.template end<Vertex>(si); ++iter)
826 {
827 Vertex* vrt = *iter;
828 UG_LOG("*\n");
829 uNew.inner_dof_indices(vrt, 0, ind);
830 DoFRef(uNew, ind[0]) = DoFRef(uOld, ind[0]);
831 }
832 };
833
834 // time step done:
835 UG_LOG("time step length: " << m_dt << "\n");
836 UG_LOG("time step no.: " << m_timestep_nr << "\n");
837 UG_LOG("time: " << m_time << "\n");
838 UG_LOG("max CFL: " << m_maxCFL << "\n");
839
840 if (m_nrOfSteps>1)
841 {
842 // Attention, uOld is overwritten
843 VecAssign(uOld,uNew);
844 };
845 };
846
847 // detach from grid
848 grid.detach_from_vertices(aScvVolume);
849 grid.detach_from_vertices(aGradient);
850
851// done
852 return true;
853}
854
858template <typename TGridFunction>
859bool FV1LevelSetDisc<TGridFunction>::calculate_vertex_vol
860(
861 TGridFunction& u,
862 aaVol& aaVolume
863)
864{
865// get domain
866 domain_type& domain = *u.domain().get();
867
868// create a FV Geometry for the dimension
869 DimFV1Geometry<dim> geo;
870
871// get position accessor
872 const position_accessor_type& aaPos = domain.position_accessor();
873
874// sum up all contributions of the sub control volumes to one vertex in an attachment
875 for(int si = 0; si < u.num_subsets(); ++si)
876 {
877 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
878 // get iterators
879 ElemIterator iter = u.template begin<ElemType>(si);
880 ElemIterator iterEnd = u.template end<ElemType>(si);
881
882 // loop elements of dimension
883 for( ;iter !=iterEnd; ++iter)
884 {
885 // get Elem
886 ElemType* elem = *iter;
887
888 size_t numVertices=elem->num_vertices();
889 // extract corner coordinates
890 for(size_t i = 0; i < numVertices; ++i)
891 coCoord[i] = aaPos[elem->vertex(i)];
892
893 // evaluate finite volume geometry
894 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
895
896
897 // loop corners
898 for (size_t i=0;i < geo.num_scv();i++)
899 {
900 // get scv for sh
901 const typename DimFV1Geometry<dim>::SCV& scv = geo.scv(i);
902
903 aaVolume[elem->vertex(i)] += scv.volume();
904 }
905 }
906 }
907
908 return true;
909}
910
916template<typename TGridFunction>
917bool FV1LevelSetDisc<TGridFunction>::compute_error(TGridFunction& numsol)
918{
919// get domain of grid function
920 domain_type& domain = *numsol.domain().get();
921
922// get grid of domain
923 grid_type& grid = *domain.grid();
924
925 // create Attachment for scv-volume size
926 ANumber aScvVolume;
927
928 // attach to grid
929 grid.attach_to_vertices(aScvVolume);
930
931 // get attachment accessor to access values
932 Grid::VertexAttachmentAccessor<ANumber> aaVolume(grid, aScvVolume);
933
934 // initialize attachment value
935 SetAttachmentValues(aaVolume, grid.vertices_begin(), grid.vertices_end(), 0);
936
937 number l1Error=0;
938 number l2Error=0;
939 number maxErr=0;
940
941 bool bRes = true;
942 // calculate scv size
943 if (! calculate_vertex_vol(numsol,aaVolume)) {UG_LOG("ERROR: gradient computation failed in compute_error function!"); };
944
945 //UG_LOG("----------------------------\n");
946
947 if(!bRes) {UG_LOG("Error while calculating CV Volume.\n"); return false;}
948 for (int si=0;si<numsol.num_subsets();++si)
949 {
950 // UG_LOG("*** " << si << "\n");
951 for(VertexConstIterator iter = numsol.template begin<Vertex>(si);
952 iter != numsol.template end<Vertex>(si); ++iter)
953 {
954 // get vertex
955 Vertex* vrt = *iter;
956 number exactVal;
957 position_accessor_type aaPos = domain.position_accessor();
958
959 // get vector holding all indices on the vertex
960 std::vector<DoFIndex> ind;
961
962 const size_t numInd = numsol.dof_indices(vrt, 0, ind);
963
964 // check indices
965 if(numInd != 1) {UG_LOG("ERROR: Wrong number of indices!"); return false;}
966
967 (*m_imDirichlet)(&exactVal,&aaPos[vrt],m_time,si,1);
968 number differ = std::abs(DoFRef(numsol, ind[0])-exactVal);
969
970 l1Error += aaVolume[vrt] * differ;
971 l2Error += aaVolume[vrt] * differ*differ;
972
973 if (m_print)
974 {
975 // if (differ>0)
976 UG_LOG("coord=" << aaPos[vrt] << " value=" << DoFRef(numsol, ind[0]) << " exact=" << exactVal << " error=" << differ << "\n");
977 };
978
979 if (differ > maxErr) maxErr = differ;
980
981 }
982 };
983 l2Error = sqrt(l2Error);
984 UG_LOG("timestep " << m_timestep_nr << " time " << m_time << "\n");
985 UG_LOG("l1 error: " << l1Error << "\n");
986 UG_LOG("l2 error: " << l2Error << "\n");
987 UG_LOG("maximum error: " << maxErr << "\n");
988 return true;
989};
990
991// next functions used for extrapolation equations as described in
992// T.D. Aslam - A partial differential equation approach to multidimensional extrapolation JCP 193 2003
993
994// compute gradient in vertices and volume of control volume in region given by sign of level set function
995template<typename TGridFunction>
996bool FV1LevelSetDisc<TGridFunction>::
997calculate_vertex_grad_vol_sign(TGridFunction& u, aaGrad& aaGradient,aaVol& aaVolume,TGridFunction& phi,int sign)
998{
999 // get domain
1000 domain_type& domain = *u.domain().get();
1001
1002 // get grid of domain
1003 typename domain_type::grid_type& grid = *domain.grid();
1004
1005 // create Multiindex
1006 std::vector<DoFIndex> multInd;
1007
1008 // create a FV Geometry for the dimension
1009 DimFV1Geometry<dim> geo;
1010
1011 // hard code function (fct=0)
1012 //\todo: generalize
1013 size_t fct=0;
1014
1015 // initialize attachment value
1016 SetAttachmentValues(aaVolume, grid.vertices_begin(), grid.vertices_end(), 0);
1017 SetAttachmentValues(aaGradient, grid.vertices_begin(), grid.vertices_end(), 0);
1018
1019 // coord and vertex array
1020 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
1021 Vertex* vVrt[domain_traits<dim>::MaxNumVerticesOfElem];
1022
1023 // sum up all contributions of the sub control volumes to one vertex in an attachment
1024 for(int si = 0; si < domain.subset_handler()->num_subsets(); ++si)
1025 {
1026 if (m_dirichlet_sg.size()!=0) if (m_dirichlet_sg.contains(si)) continue;
1027 if (m_neumann_sg.size()!=0) if (m_neumann_sg.contains(si)) continue;
1028 if (m_inactive_sg.size()!=0) if (m_inactive_sg.contains(si)) continue;
1029 // get iterators
1030 ElemIterator iter = u.template begin<ElemType>(si);
1031 ElemIterator iterEnd = u.template end<ElemType>(si);
1032
1033 // loop elements of dimension
1034 for( ;iter !=iterEnd; ++iter)
1035 {
1036 // get Elem
1037 ElemType* elem = *iter;
1038
1039 // get position accessor
1040 const position_accessor_type& aaPos = domain.position_accessor();
1041
1042 // get vertices and extract corner coordinates
1043 const size_t numVertices = elem->num_vertices();
1044 for(size_t i = 0; i < numVertices; ++i)
1045 {
1046 vVrt[i] = elem->vertex(i);
1047 coCoord[i] = aaPos[vVrt[i]];
1048 };
1049
1050 // evaluate finite volume geometry
1051 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
1052
1053 //UG_LOG("Num Verts loaded: "<<vVrt.size()<<"\n");
1054 //UG_LOG("Num SCV computed: "<<geo.num_scv()<<"\n");
1055
1056 number uValue[domain_traits<dim>::MaxNumVerticesOfElem];
1057
1058 size_t noc = geo.num_scv();
1059 bool rightsign=true;
1060 for (size_t i=0;i < noc;i++)
1061 {
1062 // get indices of function fct on vertex
1063 u.inner_dof_indices(vVrt[i], fct, multInd);
1064
1065 // read value of index from vector
1066 uValue[i]=DoFRef(u, multInd[0]);
1067 if (sign==-1)
1068 {
1069 if (DoFRef(phi, multInd[0])>0)
1070 {
1071 rightsign = false;
1072 break;
1073 }
1074 };
1075 if (sign==1)
1076 {
1077 if (DoFRef(phi, multInd[0])>0)
1078 {
1079 rightsign = false;
1080 break;
1081 };
1082 };
1083 // debug log
1084 //UG_LOG("corner " << i << " " << uValue[i] << "\n");
1085 }
1086 if (! rightsign) continue;
1087 for (size_t i=0;i < noc;i++)
1088 {
1089 // get indices of function fct on vertex
1090 u.inner_dof_indices(vVrt[i], fct, multInd);
1091
1092 // read value of index from vector
1093 uValue[i]=DoFRef(u, multInd[0]);
1094
1095 // debug log
1096 //UG_LOG("corner " << i << " " << uValue[i] << "\n");
1097 }
1098
1099 // storage for global gradient
1100 MathVector<dim> globalGrad;
1101
1102 // loop corners
1103 for (size_t i=0;i < noc;i++)
1104 {
1105 // get scv for sh
1106 const typename DimFV1Geometry<dim>::SCV& scv = geo.scv(i);
1107
1108 // debug log
1109 //UG_LOG("gradient for corner " << i << "\n");
1110
1111 // reset global gradient
1112 globalGrad = 0.0;
1113
1114 // sum up gradients of shape functions in corner
1115 for(size_t sh = 0 ; sh < noc; ++sh)
1116 {
1117 //UG_LOG("local grad " << sh << " : " << scv.local_grad(sh) << "\n");
1118 //UG_LOG("unscaled global grad " << sh << " = " << scv.global_grad(sh) << "\n");
1119 //UG_LOG("uvalue(" << sh << ") =" << uValue[sh] << "\n");
1120
1121 VecScaleAppend(globalGrad, uValue[sh], scv.global_grad(sh));
1122 }
1123
1124 // volume of scv
1125 number vol = scv.volume();
1126
1127 //UG_LOG("*** global grad " << i << ": " << globalGrad << "\n");
1128
1129 // scale gradient by volume
1130 globalGrad *= vol;
1131
1132 // add both values to attachements
1133 aaGradient[vVrt[i]] += globalGrad;
1134 aaVolume[vVrt[i]] += vol;
1135 };
1136 }
1137 }
1138
1139 // int count=0;
1140 //position_accessor_type aaPos = u.domain()->position_accessor();
1141 for (int si=0;si < u.num_subsets();++si)
1142 {
1143 //UG_LOG("si " << si << "\n");
1144 for(VertexConstIterator iter = u.template begin<Vertex>(si);
1145 iter != u.template end<Vertex>(si); ++iter)
1146 {
1147 // get vertex
1148 Vertex* vrt = *iter;
1149 if (aaVolume[vrt]!=0){
1150 (aaGradient[vrt]) /= aaVolume[vrt];
1151 };
1152 //exact[0] = cos(coord[0]);// 6*coord[0];
1153 //exact[1] = -4*sin(coord[1]);//-4*coord[1];
1154 //number gError = sqrt( (exact[0]-aaGradient[vrt][0])*(exact[0]-aaGradient[vrt][0]) + (exact[1]-aaGradient[vrt][1])*(exact[1]-aaGradient[vrt][1]) );
1155 //UG_LOG(count << "[ " << coord[0] << "," << coord[1] << " ] vol= " << aaVolume[vrt] << " " << "grad= ["
1156 // << aaGradient[vrt][0] << "," << aaGradient[vrt][1] << "] exact grad = ["
1157 // << exact[0] << "," << exact[1] << "] error: " << gError << "\n");
1158 // count++;
1159 }
1160 }
1161
1162 //UG_LOG("#*#*#*#\n");
1163 return true;
1164}
1165
1166// compute normal given by \frac{\nabla \phi}{|\nabla \phi|}
1167template<typename TGridFunction>
1168bool FV1LevelSetDisc<TGridFunction>::compute_normal(TGridFunction& vx,TGridFunction& vy,TGridFunction& u)
1169{
1170 // get grid
1171 typename domain_type::grid_type& grid = *u.domain()->grid();
1172
1173 // create Attachment for scv-volume size
1174 ANumber aScvVolume;
1175
1176 // typedef of gradient attachment
1177 typedef Attachment<MathVector<dim> > AGradient;
1178
1179 // create Attachment for gradient
1180 AGradient aGradient;
1181
1182 // attach to grid
1183 grid.attach_to_vertices(aScvVolume);
1184 grid.attach_to_vertices(aGradient);
1185
1186 // get attachment accessor to access values
1187 Grid::VertexAttachmentAccessor<ANumber> aaVolume(grid, aScvVolume);
1188 Grid::VertexAttachmentAccessor<AGradient> aaGradient(grid, aGradient);
1189
1190 // initialize attachment value
1191 SetAttachmentValues(aaVolume, grid.vertices_begin(), grid.vertices_end(), 0);
1192 SetAttachmentValues(aaGradient, grid.vertices_begin(), grid.vertices_end(), 0);
1193
1194 MathVector<dim> coord;
1195 //position_accessor_type aaPos = u.domain()->position_accessor();
1196
1197 // get vector holding all indices on the vertex
1198 std::vector<DoFIndex> ind;
1199
1200 // read indices on vertex
1201 // calculate scv size and gradient
1202 if (!calculate_vertex_grad_vol(u,aaGradient, aaVolume)) {UG_LOG("ERROR: gradient computation failed!"); return false;};
1203 if (m_limiter)
1204 limit_grad(u,aaGradient);
1205 // SetAttachmentValues(aaGradient, grid.vertices_begin(), grid.vertices_end(), 0); for debug set gradient to 0
1206 for (int si=0;si<u.num_subsets();++si)
1207 {
1208 if (m_dirichlet_sg.size()!=0) if (m_dirichlet_sg.contains(si)) continue;
1209 if (m_neumann_sg.size()!=0) if (m_neumann_sg.contains(si)) continue;
1210 if (m_inactive_sg.size()!=0) if (m_inactive_sg.contains(si)) continue;
1211 VertexConstIterator iter = u.template begin<Vertex>(si);
1212 VertexConstIterator iterEnd = u.template end<Vertex>(si);
1213 for (;iter != iterEnd; ++iter)
1214 {
1215 Vertex* vrt = *iter;
1216 u.inner_dof_indices(vrt, 0, ind);
1217 number vnorm = VecLength(aaGradient[vrt]);
1218 if (vnorm>1e-15)
1219 {
1220 DoFRef(vx, ind[0]) = aaGradient[vrt][0]/vnorm;
1221 DoFRef(vy, ind[0]) = aaGradient[vrt][1]/vnorm;
1222 } else {
1223 DoFRef(vx, ind[0]) = 0;
1224 DoFRef(vy, ind[0]) = 0;
1225 };
1226 };
1227 };
1228 return true;
1229}
1230
1231// compute directional derivative in normal direction given by \normal \cdot \nabla u (see Aslam p. 2)
1232template<typename TGridFunction>
1233bool FV1LevelSetDisc<TGridFunction>::compute_dnormal(TGridFunction& dnormal,TGridFunction& vx,TGridFunction& vy,TGridFunction& phi,TGridFunction& u)
1234{
1235 // get grid
1236 typename domain_type::grid_type& grid = *u.domain()->grid();
1237
1238 // create Attachment for scv-volume size
1239 ANumber aScvVolume;
1240
1241 // typedef of gradient attachment
1242 typedef Attachment<MathVector<dim> > AGradient;
1243
1244 // create Attachment for gradient
1245 AGradient aGradient;
1246
1247 // attach to grid
1248 grid.attach_to_vertices(aScvVolume);
1249 grid.attach_to_vertices(aGradient);
1250
1251 // get attachment accessor to access values
1252 Grid::VertexAttachmentAccessor<ANumber> aaVolume(grid, aScvVolume);
1253 Grid::VertexAttachmentAccessor<AGradient> aaGradient(grid, aGradient);
1254
1255 // initialize attachment value
1256 SetAttachmentValues(aaVolume, grid.vertices_begin(), grid.vertices_end(), 0);
1257 SetAttachmentValues(aaGradient, grid.vertices_begin(), grid.vertices_end(), 0);
1258
1259 MathVector<dim> coord;
1260
1261 // get vector holding all indices on the vertex
1262 std::vector<DoFIndex> ind;
1263 // read indices on vertex
1264 position_accessor_type aaPos = u.domain()->position_accessor();
1265//UG_LOG("-------------\n");
1266 // calculate scv size and gradient of u
1267 if (! calculate_vertex_grad_vol_sign(u,aaGradient, aaVolume,phi,-1)) {UG_LOG("ERROR: gradient computation failed!"); return false;};
1268// if (m_limiter)
1269// limit_grad(u,aaGradient);
1270 // if (calculate_vertex_grad_vol(u,aaGradient, aaVolume)==false){UG_LOG("ERROR: gradient computation failed!"); return false;};
1271 // calculate normal of phi
1272 compute_normal(vx,vy,phi);
1273 for (int si=0;si<u.num_subsets();++si)
1274 {
1275 VertexConstIterator iter = u.template begin<Vertex>(si);
1276 VertexConstIterator iterEnd = u.template end<Vertex>(si);
1277// UG_LOG("START INDEX" << si << "\n");
1278 for (;iter != iterEnd; ++iter)
1279 {
1280 Vertex* vrt = *iter;
1281 u.inner_dof_indices(vrt, 0, ind);
1282 coord = aaPos[vrt];
1283 DoFRef(dnormal, ind[0]) = DoFRef(vx, ind[0]) * aaGradient[vrt][0] + DoFRef(vy, ind[0]) * aaGradient[vrt][1];
1284// if (DoFRef(phi, ind[0]) < 0)
1285// UG_LOG("coord=(" << coord[0] << "," << coord[1] << ") exact=" << coord[0]/sqrt(coord[0]*coord[0]+coord[1]*coord[1]) << " dnormal=" << BlockRef(dnormal[ind[0][0]],ind[0][1]) << "\n" << " v=(" << BlockRef(vx[ind[0][0]],ind[0][1]) << "," << BlockRef(vy[ind[0][0]],ind[0][1]) << ") " << " grad=(" << aaGradient[vrt][0] << "," << aaGradient[vrt][1] << ")" << "\n");
1286 }
1287// UG_LOG("#\n");
1288 };
1289 return true;
1290};
1291
1292// compute directional derivative in normal direction of directional derivative in normal direction given by \normal \cdot \nabla (\normal \cdot \nabla u) (see Aslam p. 3)
1293template<typename TGridFunction>
1294bool FV1LevelSetDisc<TGridFunction>::compute_ddnormal(TGridFunction& ddnormal,TGridFunction& dnormal,TGridFunction& vx,TGridFunction& vy,TGridFunction& phi,TGridFunction& u)
1295{
1296 // get grid
1297 typename domain_type::grid_type& grid = *u.domain()->grid();
1298
1299 // create Attachment for scv-volume size
1300 ANumber aScvVolume;
1301
1302 // typedef of gradient attachment
1303 typedef Attachment<MathVector<dim> > AGradient;
1304
1305 // create Attachment for gradient
1306 AGradient aGradient;
1307
1308 // attach to grid
1309 grid.attach_to_vertices(aScvVolume);
1310 grid.attach_to_vertices(aGradient);
1311
1312 // get attachment accessor to access values
1313 Grid::VertexAttachmentAccessor<ANumber> aaVolume(grid, aScvVolume);
1314 Grid::VertexAttachmentAccessor<AGradient> aaGradient(grid, aGradient);
1315
1316 // initialize attachment value
1317 SetAttachmentValues(aaVolume, grid.vertices_begin(), grid.vertices_end(), 0);
1318 SetAttachmentValues(aaGradient, grid.vertices_begin(), grid.vertices_end(), 0);
1319
1320 MathVector<dim> coord;
1321
1322 // read indices on vertex
1323 position_accessor_type aaPos = u.domain()->position_accessor();
1324
1325 // calculate scv size and gradient of u
1326 compute_dnormal(dnormal,vx,vy,phi,u);
1327 if (! calculate_vertex_grad_vol_sign(dnormal,aaGradient, aaVolume,phi,-1)) {UG_LOG("ERROR: gradient computation failed!"); return false;};
1328 //if (m_limiter)
1329 // limit_grad(dnormal,aaGradient);
1330 // if (! calculate_vertex_grad_vol(u,aaGradient, aaVolume)) {UG_LOG("ERROR: gradient computation failed!"); return false;};
1331 // calculate normal of phi
1332
1333 std::vector<DoFIndex> ind;
1334
1335 for (int si=0;si<u.num_subsets();++si)
1336 {
1337 VertexConstIterator iter = u.template begin<Vertex>(si);
1338 VertexConstIterator iterEnd = u.template end<Vertex>(si);
1340 for (;iter != iterEnd; ++iter)
1341 {
1342 Vertex* vrt = *iter;
1343 u.inner_dof_indices(vrt, 0, ind);
1344 coord = aaPos[vrt];
1345 DoFRef(ddnormal, ind[0])
1346 = DoFRef(vx, ind[0]) * aaGradient[vrt][0]
1347 + DoFRef(vy, ind[0]) * aaGradient[vrt][1];
1348 // if (DoFRef(phi, ind[0])<0)
1349 // UG_LOG("coord=(" << coord[0] << "," << coord[1] << ") exact=" << coord[0]/sqrt(coord[0]*coord[0]+coord[1]*coord[1]) << " ddnormal=" << BlockRef(dnormal[ind[0][0]],ind[0][1]) << "\n" << " v=(" << BlockRef(vx[ind[0][0]],ind[0][1]) << "," << BlockRef(vy[ind[0][0]],ind[0][1]) << ") " << " grad=(" << aaGradient[vrt][0] << "," << aaGradient[vrt][1] << ")" << "\n");
1350 }
1352 };
1353 return true;
1354};
1355
1356// in region given by level set sign:
1357// overwrite values in parameter unew with values of uold
1358template<typename TGridFunction>
1359bool FV1LevelSetDisc<TGridFunction>::overwrite
1360(
1361 TGridFunction& unew,
1362 TGridFunction& uold,
1363 TGridFunction& phi,
1364 int sign
1365)
1366{
1367 for (int si=0;si<unew.num_subsets();++si)
1368 {
1369 for(VertexConstIterator iter = unew.template begin<Vertex>(si);
1370 iter != unew.template end<Vertex>(si); ++iter)
1371 {
1372 Vertex* vrt = * iter;
1373
1374 // read indices on vertex
1375 std::vector<DoFIndex> ind;
1376 unew.inner_dof_indices(vrt, 0, ind);
1377 number phiValue = DoFRef(phi, ind[0]);
1378 int nodeSign;
1379 if (phiValue<0) nodeSign =-1;
1380 if (phiValue>0) nodeSign = 1;
1381 if (phiValue==0) nodeSign = 0;
1382 if (nodeSign==sign) DoFRef(unew, ind[0]) = DoFRef(uold, ind[0]);
1383 }
1384 };
1385 return true;
1386};
1387
1388// in region given by level set sign:
1389// overwrite values in parameter unew with parameter value
1390template<typename TGridFunction>
1391bool FV1LevelSetDisc<TGridFunction>::overwrite
1392(
1393 TGridFunction& unew,
1394 number value,
1395 TGridFunction& phi,
1396 int sign
1397)
1398{
1399 for (int si=0;si<unew.num_subsets();++si)
1400 {
1401 for(VertexConstIterator iter = unew.template begin<Vertex>(si);
1402 iter != unew.template end<Vertex>(si); ++iter)
1403 {
1404 Vertex* vrt = * iter;
1405
1406 // read indices on vertex
1407 std::vector<MultiIndex<2> > ind;
1408 unew.inner_dof_indices(vrt, 0, ind);
1409 number phiValue = DoFRef(phi, ind[0]);
1410 int nodeSign;
1411 if (phiValue<0) nodeSign =-1;
1412 if (phiValue>0) nodeSign = 1;
1413 if (phiValue==0) nodeSign = 0;
1414 if (nodeSign==sign) DoFRef(unew, ind[0]) = value;
1415 }
1416 };
1417 return true;
1418};
1419
1420// assign subsets as given by level set function
1421// ug subset system would have to be changed for this to be useful
1422template<typename TGridFunction>
1423bool FV1LevelSetDisc<TGridFunction>::update_ls_subsets(TGridFunction& phi)
1424{
1425 // get domain of grid function
1426 domain_type& domain = *phi.domain().get();
1427
1428 // create Multiindex
1429 std::vector<DoFIndex> ind;
1430
1431 // get element iterator type
1432 m_inactive_sg.set_subset_handler(domain.subset_handler());
1433 for (int si=0;si<domain.subset_handler()->num_subsets();++si)
1434 {
1435 UG_LOG("******************* si " << si << " **********************\n");
1436 if (m_dirichlet_sg.size()!=0) if (m_dirichlet_sg.contains(si)) continue;
1437 if (m_neumann_sg.size()!=0) if (m_neumann_sg.contains(si)) continue;
1438 ElemIterator iter = phi.template begin<ElemType>(si);
1439 ElemIterator iterEnd = phi.template end<ElemType>(si);
1440 // loop elements of dimension
1441 for( ;iter !=iterEnd; ++iter)
1442 {
1443 // get Elem
1444 ElemType* elem = *iter;
1445
1446 // coord and vertex array
1447 Vertex* vVrt[domain_traits<dim>::MaxNumVerticesOfElem];
1448
1449 // get vertices and extract corner coordinates
1450 const size_t numVertices = elem->num_vertices();
1451 for(size_t i = 0; i < numVertices; ++i)
1452 {
1453 vVrt[i] = elem->vertex(i);
1454 };
1455
1456 // resize corners
1457 std::vector<MathVector<dim> > coCoord;
1458 // compute center of mass
1459 MathVector<dim> center;
1460 std::vector<MathVector<dim> > grad;
1461 center=0;
1462 int noc=elem->num_vertices();
1463 std::vector<number> phiCo(noc);
1464 for(int i = 0; i < noc; ++i)
1465 {
1466 phi.inner_dof_indices(vVrt[i], 0, ind);
1467 phiCo[i]=DoFRef(phi, ind[0]);
1468 };
1469 int firstNonzero=-1;
1470 for (int j=0;j<noc;j++)
1471 {
1472 if (phiCo[j]!=0)
1473 {
1474 firstNonzero=j;
1475 break;
1476 };
1477 };
1478 if (firstNonzero==-1) // all element nodes are on zero ls
1479 domain.subset_handler()->assign_subset(elem,m_inside_elements_si);
1480 // add nodes on zero ls to zero ls node subsetgroup
1481 for (int i=0;i<noc;i++)
1482 {
1483 if (phiCo[i]==0)
1484 {
1485 int oldindex = domain.subset_handler()->get_subset_index(vVrt[i]);
1486 if (m_dirichlet_sg.size()!=0) if (m_dirichlet_sg.contains(oldindex)) continue;
1487 if (m_neumann_sg.size()!=0) if (m_neumann_sg.contains(oldindex)) continue;
1488 domain.subset_handler()->assign_subset(vVrt[i],m_onls_nodes_si);
1489 };
1490 };
1491 bool onls=false;
1492 for (int i=firstNonzero+1;i<noc;i++)
1493 {
1494 if (phiCo[firstNonzero]*phiCo[i]<0)
1495 {
1496 onls=true;
1497 break;
1498 }
1499 };
1500 if (! onls)
1501 {
1502 if (phiCo[firstNonzero]<0)
1503 {
1504 // UG_LOG("si before " << domain.subset_handler().get_subset_index(elem));
1505 domain.subset_handler()->assign_subset(elem,m_inside_elements_si);
1506 UG_LOG("element is inside \n");
1507 // UG_LOG("si after " << domain.subset_handler().get_subset_index(elem) << "\n");
1508 // UG_LOG("-- nr of subsets: " << phi.num_subsets() << " " << m_inside_elements_si << "\n");
1509 };
1510 if (phiCo[firstNonzero]>0)
1511 {
1512 // UG_LOG("si before " << domain.subset_handler().get_subset_index(elem));
1513 domain.subset_handler()->assign_subset(elem,m_outside_elements_si);
1514 // UG_LOG("si after " << domain.subset_handler().get_subset_index(elem) << "\n");
1515 // UG_LOG("|| nr of subsets: " << phi.num_subsets() << " " << m_outside_elements_si << "\n");
1516 };
1517 };
1518 if (onls)
1519 {
1520 //UG_LOG("si before " << domain.subset_handler().get_subset_index(elem));
1521 domain.subset_handler()->assign_subset(elem,m_onls_elements_si);
1522 //UG_LOG("si after " << domain.subset_handler().get_subset_index(elem) << "\n");
1523 //UG_LOG("// nr of subsets: " << phi.num_subsets() << " " << m_onls_elements_si << "\n");
1524 for (int i=0;i<noc;i++)
1525 {
1526 int oldindex = domain.subset_handler()->get_subset_index(vVrt[i]);
1527 if (m_dirichlet_sg.size()!=0) if (m_dirichlet_sg.contains(oldindex)) continue;
1528 if (m_neumann_sg.size()!=0) if (m_neumann_sg.contains(oldindex)) continue;
1529 if (phiCo[i]<0)
1530 {
1531 domain.subset_handler()->assign_subset(vVrt[i],m_inside_nodes_si);
1532 UG_LOG("node is inside \n");
1533 };
1534 if (phiCo[i]>0)
1535 {
1536 domain.subset_handler()->assign_subset(vVrt[i],m_outside_nodes_si);
1537 };
1538 if (phiCo[i]==0)
1539 {
1540 domain.subset_handler()->assign_subset(vVrt[i],m_onls_nodes_si);
1541 }
1542 };
1543 };
1544 };
1545 };
1546/* for(int sindex = 0; sindex < phi.num_subsets(); ++sindex)
1547 {
1548 UG_LOG("si " << sindex << "\n");
1549 // get iterators
1550 ElemIterator iter = phi.template begin<ElemType>(sindex);
1551 ElemIterator iterEnd = phi.template end<ElemType>(sindex);
1552 int count=0;
1553 for( ;iter !=iterEnd; ++iter)
1554 {
1555 ElemType* elem = *iter;
1556 ++count;
1557 };
1558 UG_LOG(count << " elements in subset\n");
1559 count=0;
1560 for(VertexConstIterator iter = phi.template begin<Vertex>(sindex);
1561 iter != phi.template end<Vertex>(sindex); ++iter)
1562 {
1563 Vertex* vrt = *iter;
1564 ++count;
1565 }
1566 UG_LOG(count << " nodes in subset\n");
1567 };*/
1568 return true;
1569}
1570
1571template<typename TGridFunction>
1572bool FV1LevelSetDisc<TGridFunction>::init_ls_subsets(TGridFunction& phi)
1573{
1574 create_ls_subsets(phi);
1575 if (! update_ls_subsets(phi)) return false;
1576 return true;
1577}
1578
1579template<typename TGridFunction>
1580void FV1LevelSetDisc<TGridFunction>::create_ls_subsets(TGridFunction& phi)
1581{
1582 // get domain
1583 domain_type& domain = *phi.domain().get();
1584 UG_LOG("nr of subsets: " << domain.subset_handler()->num_subsets() << "\n");
1585 m_inside_elements_si = domain.subset_handler()->num_subsets();
1586 m_outside_elements_si = m_inside_elements_si + 1;
1587 m_onls_elements_si = m_inside_elements_si + 2;
1588 m_inside_nodes_si = m_inside_elements_si + 3;
1589 m_outside_nodes_si = m_inside_elements_si + 4;
1590 m_onls_nodes_si = m_inside_elements_si + 5;
1591 domain.subset_handler()->subset_required(phi.num_subsets()+5);
1592 UG_LOG("nr of subsets:" << domain.subset_handler()->num_subsets() << "\n");
1593};
1594
1595// for runtime testing, delete later
1596template<typename TGridFunction>
1597bool FV1LevelSetDisc<TGridFunction>::runtimetest(TGridFunction& uNew)
1598{
1599 // get domain of grid function
1600 domain_type& domain = *uNew.domain().get();
1601
1602 // create a FV Geometry for the dimension
1603 DimFV1Geometry<dim> geo;
1604
1605// get position accessor
1606 const position_accessor_type& aaPos = domain.position_accessor();
1607
1608// resize corners
1609// std::vector<MathVector<dim> > coCoord;
1610 MathVector<dim> coCoord[domain_traits<dim>::MaxNumVerticesOfElem];
1611
1612 for (int si=0;si<uNew.num_subsets();++si)
1613 {
1614 // get iterators
1615 ElemIterator iter = uNew.template begin<ElemType>(si);
1616 ElemIterator iterEnd = uNew.template end<ElemType>(si);
1617 // loop elements of dimension
1618 for( ;iter !=iterEnd; ++iter)
1619 {
1620 // get Elem
1621 ElemType* elem = *iter;
1622
1623 // extract corner coordinates
1624 const size_t numVertices = elem->num_vertices();
1625 for(size_t i = 0; i < numVertices; ++i)
1626 coCoord[i] = aaPos[elem->vertex(i)];
1627
1628 // evaluate finite volume geometry
1629 geo.update(elem, &(coCoord[0]), domain.subset_handler().get());
1630 };
1631 };
1632
1633 /*
1634 for (int si=0;si<uNew.num_subsets();++si)
1635 {
1636 // get iterators
1637 ElemIterator iter = uNew.template begin<ElemType>(si);
1638 ElemIterator iterEnd = uNew.template end<ElemType>(si);
1639 // loop elements of dimension
1640 for( ;iter !=iterEnd; ++iter)
1641 {
1642 // get Elem
1643 ElemType* elem = *iter;
1644 // get vertices of the Elem
1645 std::vector<Vertex*> vVrt;
1646 CollectVertices(vVrt, grid, elem);
1647 // get position accessor
1648 const position_accessor_type& aaPos = domain.position_accessor();
1649 // resize corners
1650 std::vector<MathVector<dim> > coCoord;
1651 // extract corner coordinates
1652 for(size_t i = 0; i < vVrt.size(); ++i)
1653 coCoord.push_back( aaPos[vVrt[i]] );
1654 };
1655 };
1656
1657
1658 for (int si=0;si<domain.subset_handler().num_subsets();++si)
1659 {
1660 // get iterators
1661 ElemIterator iter = uNew.template begin<ElemType>(si);
1662 ElemIterator iterEnd = uNew.template end<ElemType>(si);
1663 // loop elements of dimension
1664 for( ;iter !=iterEnd; ++iter)
1665 {
1666 // get Elem
1667 ElemType* elem = *iter;
1668 // get vertices of the Elem
1669 std::vector<Vertex*> vVrt;
1670 CollectVertices(vVrt, grid, elem);
1671 };
1672 };
1673
1674 MathVector<dim> coord;
1675 position_accessor_type aaPos = domain.position_accessor();
1676
1677 for (int si=0;si<uNew.num_subsets();++si)
1678 {
1679 VertexConstIterator iter = uNew.template begin<Vertex>(si);
1680 VertexConstIterator iterEnd = uNew.template end<Vertex>(si);
1681 for(;iter != iterEnd; ++iter)
1682 {
1683 // get vertex
1684 Vertex* vrt = *iter;
1685 coord = aaPos[vrt];
1686 }
1687 };*/
1688 return true;
1689}
1690
1691} // end namespace LevelSet
1692} // end namespace ug
1693
1694#endif /* LEVEL_SET_UTIL_IMPL_H_ */
attachment_util.h
SmartPtr::get
T * get()
ug::DimFV1Geometry::BF::shape
number shape(size_t sh) const
ug::DimFV1Geometry::BF::node_id
size_t node_id() const
ug::DimFV1Geometry::BF::normal
const MathVector< worldDim > & normal() const
ug::DimFV1Geometry::BF::global_grad_vector
const MathVector< worldDim > * global_grad_vector() const
ug::DimFV1Geometry::SCVF::normal
const MathVector< worldDim > & normal() const
ug::DimFV1Geometry::SCVF::global_ip
const MathVector< worldDim > & global_ip() const
ug::DimFV1Geometry::SCVF::shape
number shape(size_t sh) const
ug::DimFV1Geometry::SCV::volume
number volume() const
ug::DimFV1Geometry::SCV::global_grad
const MathVector< worldDim > & global_grad(size_t sh) const
ug::DimFV1Geometry::num_scv
size_t num_scv() const
ug::DimFV1Geometry::scv_global_ips
const MathVector< worldDim > * scv_global_ips() const
ug::DimFV1Geometry::scvf
const SCVF & scvf(size_t i) const
ug::DimFV1Geometry::num_sh
size_t num_sh() const
ug::DimFV1Geometry::scv_local_ips
const MathVector< dim > * scv_local_ips() const
ug::DimFV1Geometry::num_bf
size_t num_bf() const
ug::DimFV1Geometry::update
void update(GridObject *elem, const MathVector< worldDim > *vCornerCoords, const ISubsetHandler *ish=NULL)
ug::DimFV1Geometry::bf
const std::vector< BF > & bf(int si) const
ug::DimFV1Geometry::add_boundary_subset
void add_boundary_subset(int subsetIndex)
ug::DimFV1Geometry::scv
const SCV & scv(size_t i) const
ug::DimFV1Geometry::num_scvf
size_t num_scvf() const
ug::EdgeVertices::vertex
virtual Vertex * vertex(size_t index) const
ug::LevelSet::FV1LevelSetDisc::create_ls_subsets
void create_ls_subsets(TGridFunction &phi)
Definition level_set_impl.h:1580
ug::LevelSet::FV1LevelSetDisc::calculate_vertex_vol
bool calculate_vertex_vol(TGridFunction &u, aaVol &aaVolVolume)
fills the scvVolume attachment for all element types
Definition level_set_impl.h:860
ug::LevelSet::FV1LevelSetDisc::compute_ddnormal
bool compute_ddnormal(TGridFunction &ddnormal, TGridFunction &dnormal, TGridFunction &vx, TGridFunction &vy, TGridFunction &phi, TGridFunction &u)
Definition level_set_impl.h:1294
ug::LevelSet::FV1LevelSetDisc::update_ls_subsets
bool update_ls_subsets(TGridFunction &phi)
Definition level_set_impl.h:1423
ug::LevelSet::FV1LevelSetDisc::compute_error
bool compute_error(TGridFunction &numsol)
computes error w.r.t. the analytical solution (taken from the dirichlet bc)
Definition level_set_impl.h:917
ug::LevelSet::FV1LevelSetDisc::VertexConstIterator
TGridFunction::template traits< Vertex >::const_iterator VertexConstIterator
vertex base iterator
Definition level_set.h:89
ug::LevelSet::FV1LevelSetDisc::ElemType
TGridFunction::template dim_traits< dim >::grid_base_object ElemType
type of base grid object
Definition level_set.h:83
ug::LevelSet::FV1LevelSetDisc::assemble_element
bool assemble_element(TElem &elem, DimFV1Geometry< dim > &geo, grid_type &grid, TGridFunction &uNew, const TGridFunction &uOld, aaGrad &aaGradient, aaVol &aaVolume)
Definition level_set_impl.h:276
ug::LevelSet::FV1LevelSetDisc::grid_type
domain_type::grid_type grid_type
grid type
Definition level_set.h:71
ug::LevelSet::FV1LevelSetDisc::compute_normal
bool compute_normal(TGridFunction &vx, TGridFunction &vy, TGridFunction &u)
Definition level_set_impl.h:1168
ug::LevelSet::FV1LevelSetDisc::limit_grad
bool limit_grad(TGridFunction &uOld, aaGrad &aaGradient)
Definition level_set_impl.h:53
ug::LevelSet::FV1LevelSetDisc::runtimetest
bool runtimetest(TGridFunction &u)
Definition level_set_impl.h:1597
ug::LevelSet::FV1LevelSetDisc::position_accessor_type
domain_type::position_accessor_type position_accessor_type
type of the position accessor
Definition level_set.h:74
ug::LevelSet::FV1LevelSetDisc::ElemIterator
TGridFunction::template dim_traits< dim >::const_iterator ElemIterator
grid element iterator
Definition level_set.h:92
ug::LevelSet::FV1LevelSetDisc::assign_dirichlet
bool assign_dirichlet(TGridFunction &)
Definition level_set_impl.h:640
ug::LevelSet::FV1LevelSetDisc::compute_dnormal
bool compute_dnormal(TGridFunction &dnormal, TGridFunction &vx, TGridFunction &vy, TGridFunction &phi, TGridFunction &u)
Definition level_set_impl.h:1233
ug::LevelSet::FV1LevelSetDisc::calculate_vertex_grad_vol_sign
bool calculate_vertex_grad_vol_sign(TGridFunction &, aaGrad &, aaVol &, TGridFunction &, int)
Definition level_set_impl.h:997
ug::LevelSet::FV1LevelSetDisc::domain_type
TGridFunction::domain_type domain_type
domain type
Definition level_set.h:59
ug::LevelSet::FV1LevelSetDisc::overwrite
bool overwrite(TGridFunction &, TGridFunction &, TGridFunction &, int)
Definition level_set_impl.h:1360
ug::LevelSet::FV1LevelSetDisc::EdgeConstIterator
TGridFunction::template traits< Edge >::const_iterator EdgeConstIterator
edge iterator
Definition level_set.h:86
ug::LevelSet::FV1LevelSetDisc::init_ls_subsets
bool init_ls_subsets(TGridFunction &phi)
subset handling methods:
Definition level_set_impl.h:1572
ug::LevelSet::FV1LevelSetDisc::calculate_vertex_grad_vol
bool calculate_vertex_grad_vol(TGridFunction &u, aaGrad &aaGradient, aaVol &aaVolume)
Definition level_set_impl.h:485
ug::LevelSet::FV1LevelSetDisc::advect_lsf
bool advect_lsf(TGridFunction &uNew, TGridFunction &u)
computes the time steps of the discretization of the level-set equation
Definition level_set_impl.h:696
ug::LocalVector::resize
void resize(const LocalIndices &ind)
fv1_geom.h
ug::VecAssign
void VecAssign(GPUVector< T > &dest, const GPUVector< T > &v1)
grid
SmartPtr< TGrid > grid()
num_subsets
size_t num_subsets() const
ug::SetAttachmentValues
void SetAttachmentValues(TAttachmentAccessor &aaVal, TIter elemsBegin, TIter elemsEnd, const TVal &val)
to
bool to() const
UG_LOG
#define UG_LOG(msg)
number
double number
ug::VecLength
vector_t::value_type VecLength(const vector_t &v)
ug::VecScaleAppend
void VecScaleAppend(vector_t &vOut, typename vector_t::value_type s1, const vector_t &v1)
ug::VecAppend
void VecAppend(vector_t &vOut, const vector_t &v1)
ug::VecSubtract
void VecSubtract(vector_t &vOut, const vector_t &v, typename vector_t::value_type s)
local_finite_element_provider.h
ug::DoFRef
const number & DoFRef(const TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
ug::GetLocalVector
void GetLocalVector(LocalVector &lvec, const TVector &vec)
ug::VecScaleAdd
void VecScaleAdd(double &dest, double alpha1, const double &v1, double alpha2, const double &v2)
reference_element.h
ug::domain_traits