Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 :
34 : #include "common/util/provider.h"
35 : #include "fvho_geom.h"
36 : #include "lib_disc/reference_element/reference_element.h"
37 : #include "lib_disc/quadrature/quadrature_provider.h"
38 : #include "lib_algebra/common/operations_vec.h"
39 : #include <math.h> /* pow */
40 :
41 : namespace ug{
42 :
43 :
44 : /**
45 : * \tparam dim dimension of coordinates
46 : * \tparam TRefElem Reference element type
47 : * \tparam maxMid Maximum number of elements for all dimensions
48 : */
49 : template <int dim, typename TRefElem, int maxMid>
50 0 : static void ComputeMidPoints(const TRefElem& rRefElem,
51 : const MathVector<dim> vCorner[],
52 : MathVector<dim> vvMid[][maxMid])
53 : {
54 : // compute local midpoints for all geometric objects with 0 < d <= dim
55 0 : for(int d = 1; d <= dim; ++d)
56 : {
57 : // loop geometric objects of dimension d
58 0 : for(size_t i = 0; i < rRefElem.num(d); ++i)
59 : {
60 : // set first node
61 0 : const size_t coID0 = rRefElem.id(d, i, 0, 0);
62 0 : vvMid[d][i] = vCorner[coID0];
63 :
64 : // add corner coordinates of the corners of the geometric object
65 0 : for(size_t j = 1; j < rRefElem.num(d, i, 0); ++j)
66 : {
67 0 : const size_t coID = rRefElem.id(d, i, 0, j);
68 0 : vvMid[d][i] += vCorner[coID];
69 : }
70 :
71 : // scale for correct averaging
72 0 : vvMid[d][i] *= 1./(rRefElem.num(d, i, 0));
73 : }
74 : }
75 :
76 : // for PYRAMIDS: add midpoints of imaginary faces, edges and volumes
77 : // resulting from the division into two tetrahedra alongside x==y
78 0 : if (rRefElem.roid() == ROID_PYRAMID)
79 : {
80 0 : UG_THROW("Pyramid FV geometries are currently implemented incorrectly."
81 : " Please contact Martin Stepniewski if you see this error. ");
82 :
83 : // diagonal 2->0, diagonal 0->2
84 : VecScaleAdd(vvMid[1][rRefElem.num(1)], 0.5, vCorner[2], 0.5, vCorner[0]);
85 : VecScaleAdd(vvMid[1][rRefElem.num(1)+1], 0.5, vCorner[0], 0.5, vCorner[2]);
86 :
87 : // subface 0,1,2; subface 0,2,3; face 0,4,2; face 0,2,4
88 : vvMid[2][rRefElem.num(2)] = vCorner[0];
89 : vvMid[2][rRefElem.num(2)] += vCorner[1];
90 : vvMid[2][rRefElem.num(2)] += vCorner[2];
91 : vvMid[2][rRefElem.num(2)] *= 1.0/3.0;
92 :
93 : vvMid[2][rRefElem.num(2)+1] = vCorner[0];
94 : vvMid[2][rRefElem.num(2)+1] += vCorner[2];
95 : vvMid[2][rRefElem.num(2)+1] += vCorner[3];
96 : vvMid[2][rRefElem.num(2)+1] *= 1.0/3.0;
97 :
98 : vvMid[2][rRefElem.num(2)+2] = vCorner[0];
99 : vvMid[2][rRefElem.num(2)+2] += vCorner[4];
100 : vvMid[2][rRefElem.num(2)+2] += vCorner[2];
101 : vvMid[2][rRefElem.num(2)+2] *= 1.0/3.0;
102 :
103 : vvMid[2][rRefElem.num(2)+3] = vCorner[0];
104 : vvMid[2][rRefElem.num(2)+3] += vCorner[2];
105 : vvMid[2][rRefElem.num(2)+3] += vCorner[4];
106 : vvMid[2][rRefElem.num(2)+3] *= 1.0/3.0;
107 :
108 : // subvolume 0,1,2,4; subvolume 0,2,3,4
109 :
110 : vvMid[3][rRefElem.num(3)] = vCorner[0];
111 : vvMid[3][rRefElem.num(3)] += vCorner[1];
112 : vvMid[3][rRefElem.num(3)] += vCorner[2];
113 : vvMid[3][rRefElem.num(3)] += vCorner[4];
114 : vvMid[3][rRefElem.num(3)] *= 0.25;
115 :
116 : vvMid[3][rRefElem.num(3)+1] = vCorner[0];
117 : vvMid[3][rRefElem.num(3)+1] += vCorner[2];
118 : vvMid[3][rRefElem.num(3)+1] += vCorner[3];
119 : vvMid[3][rRefElem.num(3)+1] += vCorner[4];
120 : vvMid[3][rRefElem.num(3)+1] *= 0.25;
121 : }
122 0 : }
123 :
124 : /**
125 : * \param[in] i indicates that scvf corresponds to i'th edge of ref elem
126 : */
127 : template <typename TRefElem>
128 0 : static void ComputeSCVFMidID(const TRefElem& rRefElem,
129 : MidID vMidID[], int i)
130 : {
131 : static const int dim = TRefElem::dim;
132 :
133 0 : if (rRefElem.roid() != ROID_PYRAMID)
134 : {
135 : // set mid ids
136 : {
137 : // start at edge midpoint
138 0 : vMidID[0] = MidID(1,i);
139 :
140 : // loop up dimension
141 : if(dim == 2)
142 : {
143 0 : vMidID[1] = MidID(dim, 0); // center of element
144 : }
145 : else if (dim == 3)
146 : {
147 0 : vMidID[1] = MidID(2, rRefElem.id(1, i, 2, 0)); // side 0
148 0 : vMidID[2] = MidID(dim, 0); // center of element
149 0 : vMidID[3] = MidID(2, rRefElem.id(1, i, 2, 1)); // side 1
150 : }
151 : }
152 : }
153 : // pyramid here
154 : else
155 : {
156 0 : switch (i)
157 : {
158 : // scvf of edge 0
159 0 : case 0: vMidID[0] = MidID(1,0); // edge 0
160 0 : vMidID[1] = MidID(2,5); // subface 0/0
161 0 : vMidID[2] = MidID(3,1); // subvolume 0/0
162 0 : vMidID[3] = MidID(2,1); // face 1
163 0 : break;
164 : // scvf of edge 1
165 0 : case 1: vMidID[0] = MidID(1,1); // edge 1
166 0 : vMidID[1] = MidID(2,5); // subface 0/0
167 0 : vMidID[2] = MidID(3,1); // subvolume 0/0
168 0 : vMidID[3] = MidID(2,2); // face 2
169 0 : break;
170 : // scvf of diagonal 2->0
171 0 : case 2: vMidID[0] = MidID(1,8); // diagonal 2->0
172 0 : vMidID[1] = MidID(2,5); // subface 0/0
173 0 : vMidID[2] = MidID(3,1); // subvolume 0/0
174 0 : vMidID[3] = MidID(2,7); // face 0,4,2
175 0 : break;
176 : // scvf of edge 4 in subvolume 0/0
177 0 : case 3: vMidID[0] = MidID(1,4); // edge 4
178 0 : vMidID[1] = MidID(2,1); // face 1
179 0 : vMidID[2] = MidID(3,1); // subvolume 0/0
180 0 : vMidID[3] = MidID(2,7); // face 0,4,2
181 0 : break;
182 : // scvf of edge 5
183 0 : case 4: vMidID[0] = MidID(1,5); // edge 5
184 0 : vMidID[1] = MidID(2,2); // face 2
185 0 : vMidID[2] = MidID(3,1); // subvolume 0/0
186 0 : vMidID[3] = MidID(2,1); // face 1
187 0 : break;
188 : // scvf of edge 6 in subvolume 0/0
189 0 : case 5: vMidID[0] = MidID(1,6); // edge 6
190 0 : vMidID[1] = MidID(2,7); // face 0,4,2
191 0 : vMidID[2] = MidID(3,1); // subvolume 0/0
192 0 : vMidID[3] = MidID(2,2); // face 2
193 0 : break;
194 : // edge 0->2
195 0 : case 6: vMidID[0] = MidID(1,9); // edge 0->2
196 0 : vMidID[1] = MidID(2,6); // subface 1/0
197 0 : vMidID[2] = MidID(3,2); // subvolume 1/0
198 0 : vMidID[3] = MidID(2,8); // face 0,2,4
199 0 : break;
200 : // scvf of edge 2
201 0 : case 7: vMidID[0] = MidID(1,2); // edge 2
202 0 : vMidID[1] = MidID(2,6); // subface 1/0
203 0 : vMidID[2] = MidID(3,2); // subvolume 1/0
204 0 : vMidID[3] = MidID(2,3); // face 3
205 0 : break;
206 : // scvf of edge 3
207 0 : case 8: vMidID[0] = MidID(1,3); // edge 3
208 0 : vMidID[1] = MidID(2,6); // subface 1/0
209 0 : vMidID[2] = MidID(3,2); // subvolume 1/0
210 0 : vMidID[3] = MidID(2,4); // face 4
211 0 : break;
212 : // scvf of edge 4 in subvolume 1/0
213 0 : case 9: vMidID[0] = MidID(1,4); // edge 4
214 0 : vMidID[1] = MidID(2,8); // face 0,2,4
215 0 : vMidID[2] = MidID(3,2); // subvolume 1/0
216 0 : vMidID[3] = MidID(2,4); // face 4
217 0 : break;
218 : // scvf of edge 6 in subvolume 1/0
219 0 : case 10:vMidID[0] = MidID(1,6); // edge 6
220 0 : vMidID[1] = MidID(2,3); // face 3
221 0 : vMidID[2] = MidID(3,2); // subvolume 1/0
222 0 : vMidID[1] = MidID(2,8); // face 0,2,4
223 0 : break;
224 : // scvf of edge 7
225 0 : case 11:vMidID[0] = MidID(1,7); // edge 7
226 0 : vMidID[3] = MidID(2,4); // face 4
227 0 : vMidID[2] = MidID(3,2); // subvolume 1/0
228 0 : vMidID[1] = MidID(2,3); // face 3
229 0 : break;
230 0 : default:UG_THROW("Pyramid only has 12 SCVFs (no. 0-11), but requested no. " << i << ".");
231 : break;
232 : }
233 : }
234 0 : }
235 :
236 : /**
237 : * \param[in] i indicates that scvf corresponds to i'th corner of ref elem
238 : */
239 : template <typename TRefElem>
240 0 : static void ComputeSCVMidID(const TRefElem& rRefElem,
241 : MidID vMidID[], int i)
242 : {
243 : static const int dim = TRefElem::dim;
244 :
245 0 : if (rRefElem.roid() != ROID_PYRAMID)
246 : {
247 : if(dim == 1)
248 : {
249 0 : vMidID[0] = MidID(0, i); // set node as corner of scv
250 0 : vMidID[1] = MidID(dim, 0); // center of element
251 : }
252 : else if(dim == 2)
253 : {
254 0 : vMidID[0] = MidID(0, i); // set node as corner of scv
255 0 : vMidID[1] = MidID(1, rRefElem.id(0, i, 1, 0)); // edge 1
256 0 : vMidID[2] = MidID(dim, 0); // center of element
257 0 : vMidID[3] = MidID(1, rRefElem.id(0, i, 1, 1)); // edge 2
258 : }
259 : else if(dim == 3)
260 : {
261 0 : vMidID[0] = MidID(0, i); // set node as corner of scv
262 0 : vMidID[1] = MidID(1, rRefElem.id(0, i, 1, 1)); // edge 1
263 0 : vMidID[2] = MidID(2, rRefElem.id(0, i, 2, 0)); // face 0
264 0 : vMidID[3] = MidID(1, rRefElem.id(0, i, 1, 0)); // edge 0
265 0 : vMidID[4] = MidID(1, rRefElem.id(0, i, 1, 2)); // edge 2
266 0 : vMidID[5] = MidID(2, rRefElem.id(0, i, 2, 2)); // face 2
267 0 : vMidID[6] = MidID(dim, 0); // center of element
268 0 : vMidID[7] = MidID(2, rRefElem.id(0, i, 2, 1)); // face 1
269 : }
270 : else {UG_THROW("Dimension higher that 3 not implemented.");}
271 : }
272 : // pyramid here
273 : else
274 : {
275 0 : switch (i)
276 : {
277 : // scv of corner 0 in subvolume 0/0
278 0 : case 0: vMidID[0] = MidID(0,0); // corner 0
279 0 : vMidID[1] = MidID(1,0); // edge 0
280 0 : vMidID[2] = MidID(2,5); // subface 0/0
281 0 : vMidID[3] = MidID(1,8); // edge 2->0
282 0 : vMidID[4] = MidID(1,4); // edge 4
283 0 : vMidID[5] = MidID(2,1); // face 1
284 0 : vMidID[6] = MidID(3,1); // subvolume 0/0
285 0 : vMidID[7] = MidID(2,7); // face 0,4,2
286 0 : break;
287 : // scv of corner 1
288 0 : case 1: vMidID[0] = MidID(0,1); // corner 1
289 0 : vMidID[1] = MidID(1,1); // edge 1
290 0 : vMidID[2] = MidID(2,5); // subface 0/0
291 0 : vMidID[3] = MidID(1,0); // edge 0
292 0 : vMidID[4] = MidID(1,5); // edge 5
293 0 : vMidID[5] = MidID(2,2); // face 2
294 0 : vMidID[6] = MidID(3,1); // subvolume 0/0
295 0 : vMidID[7] = MidID(2,1); // face 1
296 0 : break;
297 : // scv of corner 2 in subvolume 0/0
298 0 : case 2: vMidID[0] = MidID(0,2); // corner 2
299 0 : vMidID[1] = MidID(1,8); // edge 2->0
300 0 : vMidID[2] = MidID(2,5); // subface 0/0
301 0 : vMidID[3] = MidID(1,1); // edge 1
302 0 : vMidID[4] = MidID(1,6); // edge 6
303 0 : vMidID[5] = MidID(2,7); // face 0,4,2
304 0 : vMidID[6] = MidID(3,1); // subvolume 0/0
305 0 : vMidID[7] = MidID(2,2); // face 2
306 0 : break;
307 : // scv of corner 4 in subvolume 0/0
308 0 : case 3: vMidID[0] = MidID(0,4); // corner 4
309 0 : vMidID[1] = MidID(1,5); // edge 5
310 0 : vMidID[2] = MidID(2,1); // face 1
311 0 : vMidID[3] = MidID(1,4); // edge 4
312 0 : vMidID[4] = MidID(1,6); // edge 6
313 0 : vMidID[5] = MidID(2,2); // face 2
314 0 : vMidID[6] = MidID(3,1); // subvolume 0/0
315 0 : vMidID[7] = MidID(2,7); // face 0,4,2
316 0 : break;
317 : // scv of corner 0 in subvolume 1/0
318 0 : case 4: vMidID[0] = MidID(0,0); // corner 0
319 0 : vMidID[1] = MidID(1,9); // edge 0->2
320 0 : vMidID[2] = MidID(2,6); // subface 1/0
321 0 : vMidID[3] = MidID(1,3); // edge 3
322 0 : vMidID[4] = MidID(1,4); // edge 4
323 0 : vMidID[5] = MidID(2,8); // face 0,2,4
324 0 : vMidID[6] = MidID(3,2); // subvolume 1/0
325 0 : vMidID[7] = MidID(2,4); // face 4
326 0 : break;
327 : // scv of corner 2 in subvolume 1/0
328 0 : case 5: vMidID[0] = MidID(0,2); // corner 2
329 0 : vMidID[1] = MidID(1,2); // edge 2
330 0 : vMidID[2] = MidID(2,6); // subface 1/0
331 0 : vMidID[3] = MidID(1,9); // edge 0->2
332 0 : vMidID[4] = MidID(1,6); // edge 6
333 0 : vMidID[5] = MidID(2,3); // face 3
334 0 : vMidID[6] = MidID(3,2); // subvolume 1/0
335 0 : vMidID[7] = MidID(2,8); // face 0,2,4
336 0 : break;
337 : // scv of corner 3
338 0 : case 6: vMidID[0] = MidID(0,3); // corner 3
339 0 : vMidID[1] = MidID(1,3); // edge 3
340 0 : vMidID[2] = MidID(2,6); // subface 1/0
341 0 : vMidID[3] = MidID(1,2); // edge 2
342 0 : vMidID[4] = MidID(1,7); // edge 7
343 0 : vMidID[5] = MidID(2,4); // face 4
344 0 : vMidID[6] = MidID(3,2); // subvolume 1/0
345 0 : vMidID[7] = MidID(2,3); // face 3
346 0 : break;
347 : // scv of corner 4 in subvolume 1/0
348 0 : case 7: vMidID[0] = MidID(0,4); // corner 4
349 0 : vMidID[1] = MidID(1,6); // edge 6
350 0 : vMidID[2] = MidID(2,8); // face 0,2,4
351 0 : vMidID[3] = MidID(1,4); // edge 4
352 0 : vMidID[4] = MidID(1,7); // edge 7
353 0 : vMidID[5] = MidID(2,3); // face 3
354 0 : vMidID[6] = MidID(3,2); // subvolume 1/0
355 0 : vMidID[7] = MidID(2,4); // face 4
356 0 : break;
357 0 : default:UG_THROW("Pyramid only has 8 SCVs (no. 0-7), but requested no. " << i << ".");
358 : break;
359 : }
360 : }
361 0 : }
362 :
363 : /**
364 : * \param[in] i indicates that scvf corresponds to i'th corner of ref elem
365 : */
366 : template <typename TRefElem>
367 0 : static void ComputeBFMidID(const TRefElem& rRefElem, int side,
368 : MidID vMidID[], int co)
369 : {
370 : static const int dim = TRefElem::dim;
371 :
372 0 : if (rRefElem.roid() != ROID_PYRAMID || side != 0)
373 : {
374 : // number of corners of side
375 0 : const int coOfSide = rRefElem.num(dim-1, side, 0);
376 :
377 : // set mid ids
378 : if(dim == 2)
379 : {
380 0 : vMidID[co%2] = MidID(0, rRefElem.id(1, side, 0, co)); // corner of side
381 0 : vMidID[(co+1)%2] = MidID(1, side); // side midpoint
382 : }
383 : else if (dim == 3)
384 : {
385 0 : vMidID[0] = MidID(0, rRefElem.id(2, side, 0, co)); // corner of side
386 0 : vMidID[1] = MidID(1, rRefElem.id(2, side, 1, co)); // edge co
387 0 : vMidID[2] = MidID(2, side); // side midpoint
388 0 : vMidID[3] = MidID(1, rRefElem.id(2, side, 1, (co -1 + coOfSide)%coOfSide)); // edge co-1
389 : }
390 : }
391 : // bottom side of pyramid here
392 : else
393 : {
394 0 : switch (co)
395 : {
396 : // bf of corner 0 in subface 0/0
397 0 : case 0: vMidID[0] = MidID(0,0); // corner 0
398 0 : vMidID[1] = MidID(1,8); // edge 2->0
399 0 : vMidID[2] = MidID(2,5); // subface 0/0
400 0 : vMidID[3] = MidID(1,0); // edge 0
401 0 : break;
402 : // bf of corner 1
403 0 : case 1: vMidID[0] = MidID(0,1); // corner 1
404 0 : vMidID[1] = MidID(1,0); // edge 0
405 0 : vMidID[2] = MidID(2,5); // subface 0/0
406 0 : vMidID[3] = MidID(1,1); // edge 1
407 0 : break;
408 : // bf of corner 2 in subvolume 0/0
409 0 : case 2: vMidID[0] = MidID(0,2); // corner 2
410 0 : vMidID[1] = MidID(1,1); // edge 1
411 0 : vMidID[2] = MidID(2,5); // subface 0/0
412 0 : vMidID[3] = MidID(1,8); // edge 2->0
413 0 : break;
414 : // bf of corner 0 in subvolume 1/0
415 0 : case 3: vMidID[0] = MidID(0,0); // corner 0
416 0 : vMidID[1] = MidID(1,3); // edge 3
417 0 : vMidID[2] = MidID(2,6); // subface 1/0
418 0 : vMidID[3] = MidID(1,9); // edge 0->2
419 0 : break;
420 : // bf of corner 2 in subvolume 1/0
421 0 : case 4: vMidID[0] = MidID(0,2); // corner 2
422 0 : vMidID[1] = MidID(1,9); // edge 0->2
423 0 : vMidID[2] = MidID(2,6); // subface 1/0
424 0 : vMidID[3] = MidID(1,2); // edge 2
425 0 : break;
426 : // bf of corner 3
427 0 : case 5: vMidID[0] = MidID(0,3); // corner 3
428 0 : vMidID[1] = MidID(1,2); // edge 2
429 0 : vMidID[2] = MidID(2,6); // subface 1/0
430 0 : vMidID[3] = MidID(1,3); // edge 3
431 0 : break;
432 0 : default:UG_THROW("Pyramid only has 6 BFs on bottom side (no. 0-5), but requested no. " << co << ".");
433 : break;
434 : }
435 : }
436 0 : }
437 :
438 : template <int dim, int maxMid>
439 0 : static void CopyCornerByMidID(MathVector<dim> vCorner[],
440 : const MidID vMidID[],
441 : MathVector<dim> vvMidPos[][maxMid],
442 : const size_t numCo)
443 : {
444 0 : for(size_t i = 0; i < numCo; ++i)
445 : {
446 0 : const size_t d = vMidID[i].dim;
447 0 : const size_t id = vMidID[i].id;
448 0 : vCorner[i] = vvMidPos[d][id];
449 : }
450 0 : }
451 :
452 : template <int dim>
453 : void ComputeMultiIndicesOfSubElement(std::vector<MathVector<dim, int> >* vvMultiIndex,
454 : bool* vIsBndElem,
455 : std::vector<int>* vElemBndSide,
456 : std::vector<size_t>* vIndex,
457 : ReferenceObjectID roid,
458 : int p);
459 :
460 : template <>
461 0 : void ComputeMultiIndicesOfSubElement<1>(std::vector<MathVector<1, int> >* vvMultiIndex,
462 : bool* vIsBndElem,
463 : std::vector<int>* vElemBndSide,
464 : std::vector<size_t>* vIndex,
465 : ReferenceObjectID roid,
466 : int p)
467 : {
468 : // switch for roid
469 : size_t se = 0;
470 0 : switch(roid)
471 : {
472 : case ROID_EDGE:
473 0 : for(int i = 0; i < p; ++i)
474 : {
475 0 : vvMultiIndex[se].resize(2);
476 : vvMultiIndex[se][0] = MathVector<1,int>(i);
477 0 : vvMultiIndex[se][1] = MathVector<1,int>(i+1);
478 :
479 : // reset bnd info
480 0 : vIsBndElem[se] = false;
481 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(3, -1);
482 :
483 0 : if(i==0)
484 : {
485 0 : vIsBndElem[se] = true;
486 0 : vElemBndSide[se][0] = 0;
487 : }
488 0 : if(i==p-1)
489 : {
490 0 : vIsBndElem[se] = true;
491 0 : vElemBndSide[se][1] = 1;
492 : }
493 0 : ++se;
494 : }
495 :
496 : {
497 0 : FlexLagrangeLSFS<ReferenceEdge> set(p);
498 0 : for(size_t s = 0; s < se; ++s)
499 : {
500 0 : vIndex[s].resize(vvMultiIndex[s].size());
501 0 : for(size_t i = 0; i < vvMultiIndex[s].size(); ++i)
502 : {
503 0 : vIndex[s][i] = set.index(vvMultiIndex[s][i]);
504 : }
505 : }
506 0 : }
507 : break;
508 0 : default: UG_THROW("ReferenceElement not found.");
509 : }
510 0 : }
511 :
512 : template <>
513 0 : void ComputeMultiIndicesOfSubElement<2>(std::vector<MathVector<2, int> >* vvMultiIndex,
514 : bool* vIsBndElem,
515 : std::vector<int>* vElemBndSide,
516 : std::vector<size_t>* vIndex,
517 : ReferenceObjectID roid,
518 : int p)
519 : {
520 : // switch for roid
521 : size_t se = 0;
522 0 : switch(roid)
523 : {
524 : case ROID_TRIANGLE:
525 0 : for(int j = 0; j < p; ++j) { // y -direction
526 0 : for(int i = 0; i < p - j; ++i) { // x - direction
527 0 : vvMultiIndex[se].resize(3);
528 : vvMultiIndex[se][0] = MathVector<2,int>(i , j );
529 0 : vvMultiIndex[se][1] = MathVector<2,int>(i+1, j );
530 0 : vvMultiIndex[se][2] = MathVector<2,int>(i , j+1);
531 :
532 : // reset bnd info
533 0 : vIsBndElem[se] = false;
534 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(3, -1);
535 :
536 0 : if(i==0) // left
537 : {
538 0 : vIsBndElem[se] = true;
539 0 : vElemBndSide[se][2] = 2;
540 : }
541 0 : if(j==0) // bottom
542 : {
543 0 : vIsBndElem[se] = true;
544 0 : vElemBndSide[se][0] = 0;
545 : }
546 0 : if(i+j==p-1) // diag
547 : {
548 0 : vIsBndElem[se] = true;
549 0 : vElemBndSide[se][1] = 1;
550 : }
551 0 : ++se;
552 : }
553 : }
554 :
555 0 : for(int j = 1; j <= p; ++j) {
556 0 : for(int i = 1; i <= p - j; ++i) {
557 0 : vvMultiIndex[se].resize(3);
558 : vvMultiIndex[se][0] = MathVector<2,int>(i , j );
559 0 : vvMultiIndex[se][1] = MathVector<2,int>(i-1, j );
560 0 : vvMultiIndex[se][2] = MathVector<2,int>(i , j-1);
561 :
562 : // reset bnd info
563 : // all inner elems
564 0 : vIsBndElem[se] = false;
565 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(3, -1);
566 0 : ++se;
567 : }
568 : }
569 :
570 : {
571 0 : FlexLagrangeLSFS<ReferenceTriangle> set(p);
572 0 : for(size_t s = 0; s < se; ++s)
573 : {
574 0 : vIndex[s].resize(vvMultiIndex[s].size());
575 0 : for(size_t i = 0; i < vvMultiIndex[s].size(); ++i)
576 : {
577 0 : vIndex[s][i] = set.index(vvMultiIndex[s][i]);
578 : }
579 : }
580 0 : }
581 :
582 0 : break;
583 : case ROID_QUADRILATERAL:
584 0 : for(int j = 0; j < p; ++j) {
585 0 : for(int i = 0; i < p; ++i) {
586 0 : vvMultiIndex[se].resize(4);
587 : vvMultiIndex[se][0] = MathVector<2,int>(i , j );
588 0 : vvMultiIndex[se][1] = MathVector<2,int>(i+1, j );
589 0 : vvMultiIndex[se][2] = MathVector<2,int>(i+1, j+1);
590 : vvMultiIndex[se][3] = MathVector<2,int>(i , j+1);
591 :
592 : // reset bnd info
593 0 : vIsBndElem[se] = false;
594 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(4, -1);
595 :
596 0 : if(i==0) // left
597 : {
598 0 : vIsBndElem[se] = true;
599 0 : vElemBndSide[se][3] = 3;
600 : }
601 0 : if(i==p-1) // right
602 : {
603 0 : vIsBndElem[se] = true;
604 0 : vElemBndSide[se][1] = 1;
605 : }
606 0 : if(j==0) // bottom
607 : {
608 0 : vIsBndElem[se] = true;
609 0 : vElemBndSide[se][0] = 0;
610 : }
611 0 : if(j==p-1) // top
612 : {
613 0 : vIsBndElem[se] = true;
614 0 : vElemBndSide[se][2] = 2;
615 : }
616 0 : ++se;
617 : }
618 : }
619 :
620 : {
621 0 : FlexLagrangeLSFS<ReferenceQuadrilateral> set(p);
622 0 : for(size_t s = 0; s < se; ++s)
623 : {
624 0 : vIndex[s].resize(vvMultiIndex[s].size());
625 0 : for(size_t i = 0; i < vvMultiIndex[s].size(); ++i)
626 : {
627 0 : vIndex[s][i] = set.index(vvMultiIndex[s][i]);
628 : }
629 : }
630 0 : }
631 0 : break;
632 0 : default: UG_THROW("ReferenceElement not found.");
633 : }
634 :
635 0 : }
636 :
637 : template <>
638 0 : void ComputeMultiIndicesOfSubElement<3>(std::vector<MathVector<3, int> >* vvMultiIndex,
639 : bool* vIsBndElem,
640 : std::vector<int>* vElemBndSide,
641 : std::vector<size_t>* vIndex,
642 : ReferenceObjectID roid,
643 : int p)
644 : {
645 : // switch for roid
646 : size_t se = 0;
647 0 : switch(roid)
648 : {
649 : case ROID_TETRAHEDRON:
650 0 : for(int k = 0; k < p; ++k) {
651 0 : for(int j = 0; j < p -k; ++j) {
652 0 : for(int i = 0; i < p -k -j; ++i) {
653 0 : vvMultiIndex[se].resize(4);
654 : vvMultiIndex[se][0] = MathVector<3,int>(i , j , k);
655 0 : vvMultiIndex[se][1] = MathVector<3,int>(i+1, j , k);
656 0 : vvMultiIndex[se][2] = MathVector<3,int>(i , j+1, k);
657 0 : vvMultiIndex[se][3] = MathVector<3,int>(i , j , k+1);
658 :
659 : // reset bnd info
660 0 : vIsBndElem[se] = false;
661 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(4, -1);
662 :
663 0 : if(i==0) // left
664 : {
665 0 : vIsBndElem[se] = true;
666 0 : vElemBndSide[se][2] = 2;
667 : }
668 0 : if(j==0) // front
669 : {
670 0 : vIsBndElem[se] = true;
671 0 : vElemBndSide[se][3] = 3;
672 : }
673 0 : if(k==0) // bottom
674 : {
675 0 : vIsBndElem[se] = true;
676 0 : vElemBndSide[se][0] = 0;
677 : }
678 0 : if(i+j+k==p-1) // diag
679 : {
680 0 : vIsBndElem[se] = true;
681 0 : vElemBndSide[se][1] = 1;
682 : }
683 0 : ++se;
684 : }
685 : }
686 : }
687 : // build 4 tetrahedrons out of the remaining octogons
688 0 : for(int k = 0; k < p; ++k) {
689 0 : for(int j = 1; j < p -k; ++j) {
690 0 : for(int i = 0; i < p -k -j; ++i) {
691 0 : if(j >= 2){
692 0 : vvMultiIndex[se].resize(4);
693 0 : vvMultiIndex[se][0] = MathVector<3,int>(i+1, j-1, k);
694 0 : vvMultiIndex[se][1] = MathVector<3,int>(i+1, j-1, k+1);
695 : vvMultiIndex[se][2] = MathVector<3,int>(i , j-1, k+1);
696 0 : vvMultiIndex[se][3] = MathVector<3,int>(i+1, j-2, k+1);
697 :
698 : // reset bnd info
699 0 : vIsBndElem[se] = false;
700 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(4, -1);
701 :
702 0 : ++se;
703 : }
704 :
705 0 : vvMultiIndex[se].resize(4);
706 : vvMultiIndex[se][0] = MathVector<3,int>(i , j , k);
707 0 : vvMultiIndex[se][1] = MathVector<3,int>(i , j , k+1);
708 0 : vvMultiIndex[se][2] = MathVector<3,int>(i , j-1, k+1);
709 0 : vvMultiIndex[se][3] = MathVector<3,int>(i+1, j-1, k+1);
710 :
711 : // reset bnd info
712 0 : vIsBndElem[se] = false;
713 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(4, -1);
714 :
715 0 : if(i==0) // left
716 : {
717 0 : vIsBndElem[se] = true;
718 0 : vElemBndSide[se][0] = 2;
719 : }
720 0 : ++se;
721 :
722 0 : vvMultiIndex[se].resize(4);
723 : vvMultiIndex[se][0] = MathVector<3,int>(i , j , k);
724 0 : vvMultiIndex[se][1] = MathVector<3,int>(i+1, j , k);
725 : vvMultiIndex[se][2] = MathVector<3,int>(i+1, j-1, k+1);
726 : vvMultiIndex[se][3] = MathVector<3,int>(i+1, j-1, k);
727 :
728 : // reset bnd info
729 0 : vIsBndElem[se] = false;
730 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(4, -1);
731 :
732 0 : if(k==0) // bottom
733 : {
734 0 : vIsBndElem[se] = true;
735 0 : vElemBndSide[se][2] = 0;
736 : }
737 0 : ++se;
738 :
739 0 : vvMultiIndex[se].resize(4);
740 : vvMultiIndex[se][0] = MathVector<3,int>(i , j , k);
741 : vvMultiIndex[se][1] = MathVector<3,int>(i+1, j-1, k);
742 : vvMultiIndex[se][2] = MathVector<3,int>(i+1, j-1, k+1);
743 : vvMultiIndex[se][3] = MathVector<3,int>(i , j-1, k+1);
744 :
745 : // reset bnd info
746 0 : vIsBndElem[se] = false;
747 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(4, -1);
748 :
749 0 : if(j==1) // front
750 : {
751 0 : vIsBndElem[se] = true;
752 0 : vElemBndSide[se][1] = 3;
753 : }
754 0 : ++se;
755 :
756 0 : vvMultiIndex[se].resize(4);
757 : vvMultiIndex[se][0] = MathVector<3,int>(i , j , k);
758 : vvMultiIndex[se][1] = MathVector<3,int>(i+1, j-1, k+1);
759 : vvMultiIndex[se][2] = MathVector<3,int>(i+1, j , k);
760 : vvMultiIndex[se][3] = MathVector<3,int>(i , j , k+1);
761 :
762 : // reset bnd info
763 0 : vIsBndElem[se] = false;
764 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(4, -1);
765 :
766 0 : if(i+j+k==p-1) // diag
767 : {
768 0 : vIsBndElem[se] = true;
769 0 : vElemBndSide[se][1] = 1;
770 : }
771 0 : ++se;
772 : }
773 : }
774 : }
775 :
776 : {
777 0 : FlexLagrangeLSFS<ReferenceTetrahedron> set(p);
778 0 : for(size_t s = 0; s < se; ++s)
779 : {
780 0 : vIndex[s].resize(vvMultiIndex[s].size());
781 0 : for(size_t i = 0; i < vvMultiIndex[s].size(); ++i)
782 : {
783 0 : vIndex[s][i] = set.index(vvMultiIndex[s][i]);
784 : }
785 : }
786 0 : }
787 0 : break;
788 :
789 : case ROID_PRISM:
790 0 : for(int k = 0; k < p; ++k) {
791 0 : for(int j = 0; j < p; ++j) {
792 0 : for(int i = 0; i < p - j; ++i) {
793 0 : vvMultiIndex[se].resize(6);
794 : vvMultiIndex[se][0] = MathVector<3,int>(i , j , k);
795 0 : vvMultiIndex[se][1] = MathVector<3,int>(i+1, j , k);
796 0 : vvMultiIndex[se][2] = MathVector<3,int>(i , j+1, k);
797 0 : vvMultiIndex[se][3] = MathVector<3,int>(i , j , k+1);
798 : vvMultiIndex[se][4] = MathVector<3,int>(i+1, j , k+1);
799 : vvMultiIndex[se][5] = MathVector<3,int>(i , j+1, k+1);
800 :
801 : // reset bnd info
802 0 : vIsBndElem[se] = false;
803 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(5, -1);
804 :
805 0 : if(i==0) // left
806 : {
807 0 : vIsBndElem[se] = true;
808 0 : vElemBndSide[se][3] = 3;
809 : }
810 0 : if(j==0) // front
811 : {
812 0 : vIsBndElem[se] = true;
813 0 : vElemBndSide[se][1] = 1;
814 : }
815 0 : if(i+j==p-1) // diag
816 : {
817 0 : vIsBndElem[se] = true;
818 0 : vElemBndSide[se][2] = 2;
819 : }
820 0 : if(k==0) // bottom
821 : {
822 0 : vIsBndElem[se] = true;
823 0 : vElemBndSide[se][0] = 0;
824 : }
825 0 : if(k==p-1) // top
826 : {
827 0 : vIsBndElem[se] = true;
828 0 : vElemBndSide[se][4] = 4;
829 : }
830 0 : ++se;
831 : }
832 : }
833 : }
834 0 : for(int k = 0; k < p; ++k) {
835 0 : for(int j = 1; j <= p; ++j) {
836 0 : for(int i = 1; i <= p - j; ++i) {
837 0 : vvMultiIndex[se].resize(6);
838 : vvMultiIndex[se][0] = MathVector<3,int>(i , j ,k);
839 0 : vvMultiIndex[se][1] = MathVector<3,int>(i-1, j ,k);
840 0 : vvMultiIndex[se][2] = MathVector<3,int>(i , j-1,k);
841 0 : vvMultiIndex[se][3] = MathVector<3,int>(i , j ,k+1);
842 : vvMultiIndex[se][4] = MathVector<3,int>(i-1, j ,k+1);
843 : vvMultiIndex[se][5] = MathVector<3,int>(i , j-1,k+1);
844 :
845 : // reset bnd info
846 0 : vIsBndElem[se] = false;
847 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(5, -1);
848 :
849 0 : if(k==0) // bottom
850 : {
851 0 : vIsBndElem[se] = true;
852 0 : vElemBndSide[se][0] = 0;
853 : }
854 0 : if(k==p-1) // top
855 : {
856 0 : vIsBndElem[se] = true;
857 0 : vElemBndSide[se][4] = 4;
858 : }
859 0 : ++se;
860 : }
861 : }
862 : }
863 :
864 : {
865 0 : FlexLagrangeLSFS<ReferencePrism> set(p);
866 0 : for(size_t s = 0; s < se; ++s)
867 : {
868 0 : vIndex[s].resize(vvMultiIndex[s].size());
869 0 : for(size_t i = 0; i < vvMultiIndex[s].size(); ++i)
870 : {
871 0 : vIndex[s][i] = set.index(vvMultiIndex[s][i]);
872 : }
873 : }
874 0 : }
875 :
876 0 : break;
877 : case ROID_HEXAHEDRON:
878 0 : for(int k = 0; k < p; ++k) {
879 0 : for(int j = 0; j < p; ++j) {
880 0 : for(int i = 0; i < p; ++i) {
881 0 : vvMultiIndex[se].resize(8);
882 : vvMultiIndex[se][0] = MathVector<3,int>(i , j , k);
883 0 : vvMultiIndex[se][1] = MathVector<3,int>(i+1, j , k);
884 0 : vvMultiIndex[se][2] = MathVector<3,int>(i+1, j+1, k);
885 : vvMultiIndex[se][3] = MathVector<3,int>(i , j+1, k);
886 0 : vvMultiIndex[se][4] = MathVector<3,int>(i , j , k+1);
887 : vvMultiIndex[se][5] = MathVector<3,int>(i+1, j , k+1);
888 : vvMultiIndex[se][6] = MathVector<3,int>(i+1, j+1, k+1);
889 : vvMultiIndex[se][7] = MathVector<3,int>(i , j+1, k+1);
890 :
891 : // reset bnd info
892 0 : vIsBndElem[se] = false;
893 0 : vElemBndSide[se].clear(); vElemBndSide[se].resize(6, -1);
894 :
895 0 : if(i==0) // left
896 : {
897 0 : vIsBndElem[se] = true;
898 0 : vElemBndSide[se][4] = 4;
899 : }
900 0 : if(i==p-1) //right
901 : {
902 0 : vIsBndElem[se] = true;
903 0 : vElemBndSide[se][2] = 2;
904 : }
905 0 : if(j==0) // front
906 : {
907 0 : vIsBndElem[se] = true;
908 0 : vElemBndSide[se][1] = 1;
909 : }
910 0 : if(j==p-1) // back
911 : {
912 0 : vIsBndElem[se] = true;
913 0 : vElemBndSide[se][3] = 3;
914 : }
915 0 : if(k==0) // bottom
916 : {
917 0 : vIsBndElem[se] = true;
918 0 : vElemBndSide[se][0] = 0;
919 : }
920 0 : if(k==p-1) // top
921 : {
922 0 : vIsBndElem[se] = true;
923 0 : vElemBndSide[se][5] = 5;
924 : }
925 0 : ++se;
926 : }
927 : }
928 : }
929 :
930 : {
931 0 : FlexLagrangeLSFS<ReferenceHexahedron> set(p);
932 0 : for(size_t s = 0; s < se; ++s)
933 : {
934 0 : vIndex[s].resize(vvMultiIndex[s].size());
935 0 : for(size_t i = 0; i < vvMultiIndex[s].size(); ++i)
936 : {
937 0 : vIndex[s][i] = set.index(vvMultiIndex[s][i]);
938 : }
939 : }
940 0 : }
941 :
942 0 : break;
943 0 : default: UG_THROW("ReferenceElement not found.");
944 : }
945 :
946 0 : }
947 :
948 :
949 : ////////////////////////////////////////////////////////////////////////////////
950 : ////////////////////////////////////////////////////////////////////////////////
951 : // FV Geometry for Reference Element Type (all order, FVHO)
952 : ////////////////////////////////////////////////////////////////////////////////
953 : ////////////////////////////////////////////////////////////////////////////////
954 :
955 : template <int TOrder, typename TElem, int TWorldDim, int TQuadOrder>
956 0 : FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>::
957 : FVGeometry()
958 0 : : m_pElem(NULL), m_rRefElem(Provider<ref_elem_type>::get()),
959 0 : m_rTrialSpace(Provider<local_shape_fct_set_type>::get()),
960 0 : m_rSCVFQuadRule(Provider<scvf_quad_rule_type>::get()),
961 0 : m_rSCVQuadRule(Provider<scv_quad_rule_type>::get())
962 : {
963 0 : update_local_data();
964 0 : }
965 :
966 : template <int TOrder, typename TElem, int TWorldDim, int TQuadOrder>
967 0 : void FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>::
968 : update_local(ReferenceObjectID roid, const LFEID& lfeID,
969 : size_t quadOrder)
970 : {
971 0 : if(roid != geometry_traits<TElem>::REFERENCE_OBJECT_ID)
972 : {
973 0 : UG_THROW("FVGeometry::update: Geometry only for "
974 : <<geometry_traits<TElem>::REFERENCE_OBJECT_ID<<", but "
975 : <<roid<<" requested.");
976 : }
977 0 : if(lfeID.type() != LFEID::LAGRANGE)
978 : {
979 0 : UG_THROW("FVGeometry::update: Geometry only for shape type"
980 : <<"Lagrange"<<", but "<<lfeID.type()<<" requested.");
981 : }
982 0 : if(lfeID.order() != TOrder)
983 : {
984 0 : UG_THROW("FVGeometry::update: Geometry only for shape order"
985 : <<TOrder<<", but "<<lfeID.order()<<" requested.");
986 : }
987 0 : if(quadOrder > TQuadOrder)
988 : {
989 0 : UG_THROW("FVGeometry::update: Geometry only for scvf integration order "
990 : << TQuadOrder<<", but order "<<quadOrder<<" requested.");
991 : }
992 0 : }
993 :
994 :
995 : template <int TOrder, typename TElem, int TWorldDim, int TQuadOrder>
996 0 : void FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>::
997 : update_local_data()
998 : {
999 : // get reference object id
1000 : ReferenceObjectID roid = ref_elem_type::REFERENCE_OBJECT_ID;
1001 :
1002 : // determine corners of sub elements
1003 : bool vIsBndElem[numSubElem];
1004 0 : std::vector<int> vElemBndSide[numSubElem];
1005 0 : std::vector<size_t> vIndex[numSubElem];
1006 0 : std::vector<MathVector<dim,int> > vMultiIndex[numSubElem];
1007 :
1008 0 : ComputeMultiIndicesOfSubElement<dim>(vMultiIndex, vIsBndElem,
1009 : vElemBndSide, vIndex, roid, p);
1010 :
1011 : // directions of counting
1012 0 : MathVector<dim> direction[dim];
1013 0 : for(int i = 0; i < dim; ++i){direction[i] = 0.0; direction[i][i] = 1.0;}
1014 :
1015 0 : for(size_t se = 0; se < numSubElem; ++se)
1016 : {
1017 0 : for(int co = 0; co < ref_elem_type::numCorners; ++co)
1018 : {
1019 : // compute corners of sub elem in local coordinates
1020 : MathVector<dim> pos; pos = 0.0;
1021 0 : for(int i = 0; i < dim; ++i)
1022 : {
1023 0 : const number frac = vMultiIndex[se][co][i] / ((number)p);
1024 : VecScaleAppend(pos, frac, direction[i]);
1025 : }
1026 : m_vSubElem[se].vvLocMid[0][co] = pos;
1027 :
1028 : // get multi index for corner
1029 0 : m_vSubElem[se].vDoFID[co] = vIndex[se][co];
1030 : }
1031 :
1032 : // remember if boundary element
1033 0 : m_vSubElem[se].isBndElem = vIsBndElem[se];
1034 :
1035 : // remember boundary sides
1036 0 : m_vSubElem[se].vElemBndSide = vElemBndSide[se];
1037 : }
1038 :
1039 : // compute mid points for all sub elements
1040 0 : for(size_t se = 0; se < numSubElem; ++se)
1041 : ComputeMidPoints<dim, ref_elem_type, maxMid>
1042 0 : (m_rRefElem, m_vSubElem[se].vvLocMid[0], m_vSubElem[se].vvLocMid);
1043 :
1044 : // set up local informations for SubControlVolumeFaces (scvf)
1045 : // each scvf is associated to one edge of the sub-element
1046 0 : for(size_t i = 0; i < num_scvf(); ++i)
1047 : {
1048 : // get corresponding subelement
1049 0 : const size_t se = i / numSCVFPerSubElem;
1050 0 : const size_t locSCVF = i % numSCVFPerSubElem;
1051 :
1052 : // this scvf separates the given nodes
1053 0 : const size_t locFrom = m_rRefElem.id(1, locSCVF, 0, 0);
1054 0 : const size_t locTo = m_rRefElem.id(1, locSCVF, 0, 1);
1055 :
1056 0 : m_vSCVF[i].From = m_vSubElem[se].vDoFID[locFrom];
1057 0 : m_vSCVF[i].To = m_vSubElem[se].vDoFID[locTo];
1058 :
1059 : // compute mid ids of the scvf
1060 0 : ComputeSCVFMidID(m_rRefElem, m_vSCVF[i].vMidID, locSCVF);
1061 :
1062 : // copy local corners of scvf
1063 : CopyCornerByMidID<dim, maxMid>
1064 0 : (m_vSCVF[i].vLocPos, m_vSCVF[i].vMidID, m_vSubElem[se].vvLocMid, SCVF::numCo);
1065 :
1066 : // compute integration points
1067 0 : m_vSCVF[i].vWeight = m_rSCVFQuadRule.weights();
1068 0 : ReferenceMapping<scvf_type, dim> map(m_vSCVF[i].vLocPos);
1069 0 : for(size_t ip = 0; ip < m_rSCVFQuadRule.size(); ++ip)
1070 0 : map.local_to_global(m_vSCVF[i].vLocalIP[ip], m_rSCVFQuadRule.point(ip));
1071 : }
1072 :
1073 :
1074 : // set up local informations for SubControlVolumes (scv)
1075 : // each scv is associated to one corner of the sub-element
1076 0 : for(size_t i = 0; i < num_scv(); ++i)
1077 : {
1078 : // get corresponding subelement
1079 0 : const size_t se = i / numSCVPerSubElem;
1080 0 : const size_t locSCV = i % numSCVPerSubElem;
1081 :
1082 : // store associated node
1083 0 : m_vSCV[i].nodeId = m_vSubElem[se].vDoFID[locSCV];;
1084 :
1085 : // compute mid ids scv
1086 0 : ComputeSCVMidID(m_rRefElem, m_vSCV[i].vMidID, locSCV);
1087 :
1088 : // copy local corners of scv
1089 : CopyCornerByMidID<dim, maxMid>
1090 0 : (m_vSCV[i].vLocPos, m_vSCV[i].vMidID, m_vSubElem[se].vvLocMid, m_vSCV[i].num_corners());
1091 :
1092 : // compute integration points
1093 0 : m_vSCV[i].vWeight = m_rSCVQuadRule.weights();
1094 0 : ReferenceMapping<scv_type, dim> map(m_vSCV[i].vLocPos);
1095 0 : for(size_t ip = 0; ip < m_rSCVQuadRule.size(); ++ip)
1096 0 : map.local_to_global(m_vSCV[i].vLocalIP[ip], m_rSCVQuadRule.point(ip));
1097 : }
1098 :
1099 : /////////////////////////
1100 : // Shapes and Derivatives
1101 : /////////////////////////
1102 :
1103 0 : for(size_t i = 0; i < num_scvf(); ++i)
1104 0 : for(size_t ip = 0; ip < m_vSCVF[i].num_ip(); ++ip)
1105 : {
1106 0 : m_rTrialSpace.shapes(&(m_vSCVF[i].vvShape[ip][0]), m_vSCVF[i].local_ip(ip));
1107 0 : m_rTrialSpace.grads(&(m_vSCVF[i].vvLocalGrad[ip][0]), m_vSCVF[i].local_ip(ip));
1108 : }
1109 :
1110 0 : for(size_t i = 0; i < num_scv(); ++i)
1111 0 : for(size_t ip = 0; ip < m_vSCV[i].num_ip(); ++ip)
1112 : {
1113 0 : m_rTrialSpace.shapes(&(m_vSCV[i].vvShape[ip][0]), m_vSCV[i].local_ip(ip));
1114 0 : m_rTrialSpace.grads(&(m_vSCV[i].vvLocalGrad[ip][0]), m_vSCV[i].local_ip(ip));
1115 : }
1116 :
1117 : // copy ip positions in a list for Sub Control Volumes Faces (SCVF)
1118 : size_t allIP = 0;
1119 0 : for(size_t i = 0; i < num_scvf(); ++i)
1120 0 : for(size_t ip = 0; ip < m_vSCVF[i].num_ip(); ++ip)
1121 0 : m_vLocSCVF_IP[allIP++] = scvf(i).local_ip(ip);
1122 :
1123 : // copy ip positions in a list for Sub Control Volumes (SCV)
1124 : allIP = 0;
1125 0 : for(size_t i = 0; i < num_scv(); ++i)
1126 0 : for(size_t ip = 0; ip < m_vSCV[i].num_ip(); ++ip)
1127 0 : m_vLocSCV_IP[allIP++] = scv(i).local_ip(ip);
1128 0 : }
1129 :
1130 :
1131 : /// update data for given element
1132 : template <int TOrder, typename TElem, int TWorldDim, int TQuadOrder>
1133 0 : void FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>::
1134 : update(GridObject* elem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
1135 : {
1136 : UG_ASSERT(dynamic_cast<TElem*>(elem) != NULL, "Wrong element type.");
1137 : TElem* pElem = static_cast<TElem*>(elem);
1138 :
1139 : // If already update for this element, do nothing
1140 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
1141 :
1142 : // update reference mapping
1143 0 : m_rMapping.update(vCornerCoords);
1144 :
1145 : // compute global informations for scvf
1146 0 : for(size_t i = 0; i < num_scvf(); ++i)
1147 : {
1148 : // map local corners of scvf to global
1149 0 : for(size_t co = 0; co < m_vSCVF[i].num_corners(); ++co)
1150 0 : m_rMapping.local_to_global(m_vSCVF[i].vGloPos[co], m_vSCVF[i].vLocPos[co]);
1151 :
1152 : // map local ips of scvf to global
1153 0 : for(size_t ip = 0; ip < m_vSCVF[i].num_ip(); ++ip)
1154 0 : m_rMapping.local_to_global(m_vSCVF[i].vGlobalIP[ip], m_vSCVF[i].local_ip(ip));
1155 :
1156 : // normal on scvf
1157 0 : traits::NormalOnSCVF(m_vSCVF[i].Normal, m_vSCVF[i].vGloPos, vCornerCoords);
1158 0 : VecNormalize(m_vSCVF[i].Normal, m_vSCVF[i].Normal);
1159 :
1160 0 : ReferenceMapping<scvf_type, worldDim> map(&m_vSCVF[i].vGloPos[0]);
1161 0 : for(size_t ip = 0; ip < m_rSCVFQuadRule.size(); ++ip)
1162 0 : m_vSCVF[i].vDetJMap[ip] = map.sqrt_gram_det(m_rSCVFQuadRule.point(ip));
1163 : }
1164 :
1165 : // compute size of scv
1166 0 : for(size_t i = 0; i < num_scv(); ++i)
1167 : {
1168 : // map local corners of scvf to global
1169 0 : for(size_t co = 0; co < m_vSCV[i].num_corners(); ++co)
1170 0 : m_rMapping.local_to_global(m_vSCV[i].vGloPos[co], m_vSCV[i].vLocPos[co]);
1171 :
1172 : // map local ips of scvf to global
1173 0 : for(size_t ip = 0; ip < m_vSCV[i].num_ip(); ++ip)
1174 0 : m_rMapping.local_to_global(m_vSCV[i].vGlobalIP[ip], m_vSCV[i].local_ip(ip));
1175 : }
1176 :
1177 : // if mapping is linear, compute jacobian only once and copy
1178 : if(ReferenceMapping<ref_elem_type, worldDim>::isLinear)
1179 : {
1180 : MathMatrix<worldDim,dim> JtInv;
1181 0 : m_rMapping.jacobian_transposed_inverse(JtInv, m_vSCVF[0].local_ip(0));
1182 0 : const number detJ = m_rMapping.sqrt_gram_det(m_vSCVF[0].local_ip(0));
1183 0 : for(size_t i = 0; i < num_scvf(); ++i)
1184 0 : for(size_t ip = 0; ip < scvf(i).num_ip(); ++ip)
1185 : {
1186 : m_vSCVF[i].vJtInv[ip] = JtInv;
1187 0 : m_vSCVF[i].vDetJ[ip] = detJ;
1188 : }
1189 :
1190 0 : for(size_t i = 0; i < num_scv(); ++i)
1191 0 : for(size_t ip = 0; ip < scv(i).num_ip(); ++ip)
1192 : {
1193 : m_vSCV[i].vJtInv[ip] = JtInv;
1194 0 : m_vSCV[i].vDetJ[ip] = detJ;
1195 : }
1196 : }
1197 : // else compute jacobian for each integration point
1198 : else
1199 : {
1200 0 : for(size_t i = 0; i < num_scvf(); ++i)
1201 0 : for(size_t ip = 0; ip < m_vSCVF[i].num_ip(); ++ip)
1202 : {
1203 0 : m_rMapping.jacobian_transposed_inverse(m_vSCVF[i].vJtInv[ip], m_vSCVF[i].local_ip(ip));
1204 0 : m_vSCVF[i].vDetJ[ip] = m_rMapping.sqrt_gram_det(m_vSCVF[i].local_ip(ip));
1205 : }
1206 :
1207 0 : for(size_t i = 0; i < num_scv(); ++i)
1208 0 : for(size_t ip = 0; ip < m_vSCV[i].num_ip(); ++ip)
1209 : {
1210 0 : m_rMapping.jacobian_transposed_inverse(m_vSCV[i].vJtInv[ip], m_vSCV[i].local_ip(ip));
1211 0 : m_vSCV[i].vDetJ[ip] = m_rMapping.sqrt_gram_det(m_vSCV[i].local_ip(ip));
1212 : }
1213 : }
1214 :
1215 0 : for(size_t i = 0; i < num_scv(); ++i)
1216 : {
1217 0 : ReferenceMapping<scv_type, worldDim> map(&m_vSCV[i].vGloPos[0]);
1218 0 : for(size_t ip = 0; ip < m_rSCVQuadRule.size(); ++ip)
1219 0 : m_vSCV[i].vDetJMap[ip] = map.sqrt_gram_det(m_rSCVQuadRule.point(ip));
1220 : }
1221 :
1222 : // compute global gradients
1223 0 : for(size_t i = 0; i < num_scvf(); ++i)
1224 0 : for(size_t ip = 0; ip < scvf(i).num_ip(); ++ip)
1225 0 : for(size_t sh = 0 ; sh < nsh; ++sh)
1226 0 : MatVecMult(m_vSCVF[i].vvGlobalGrad[ip][sh], m_vSCVF[i].vJtInv[ip], m_vSCVF[i].vvLocalGrad[ip][sh]);
1227 :
1228 0 : for(size_t i = 0; i < num_scv(); ++i)
1229 0 : for(size_t ip = 0; ip < scv(i).num_ip(); ++ip)
1230 0 : for(size_t sh = 0 ; sh < nsh; ++sh)
1231 0 : MatVecMult(m_vSCV[i].vvGlobalGrad[ip][sh], m_vSCV[i].vJtInv[ip], m_vSCV[i].vvLocalGrad[ip][sh]);
1232 :
1233 : // Copy ip pos in list for SCVF
1234 : size_t allIP = 0;
1235 0 : for(size_t i = 0; i < num_scvf(); ++i)
1236 0 : for(size_t ip = 0; ip < scvf(i).num_ip(); ++ip)
1237 0 : m_vGlobSCVF_IP[allIP++] = scvf(i).global_ip(ip);
1238 :
1239 : allIP = 0;
1240 0 : for(size_t i = 0; i < num_scv(); ++i)
1241 0 : for(size_t ip = 0; ip < scv(i).num_ip(); ++ip)
1242 0 : m_vGlobSCV_IP[allIP++] = scv(i).global_ip(ip);
1243 :
1244 : // if no boundary subsets required, return
1245 0 : if(num_boundary_subsets() == 0 || ish == NULL) return;
1246 0 : else update_boundary_faces(pElem, vCornerCoords, ish);
1247 : }
1248 :
1249 : template <int TOrder, typename TElem, int TWorldDim, int TQuadOrder>
1250 0 : void FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>::
1251 : update_boundary_faces(GridObject* elem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
1252 : {
1253 : UG_ASSERT(dynamic_cast<TElem*>(elem) != NULL, "Wrong element type.");
1254 : TElem* pElem = static_cast<TElem*>(elem);
1255 :
1256 : // get grid
1257 0 : Grid& grid = *(ish->grid());
1258 :
1259 : // vector of subset indices of side
1260 : std::vector<int> vSubsetIndex;
1261 :
1262 : // get subset indices for sides (i.e. edge in 2d, faces in 3d)
1263 : if(dim == 1) {
1264 : std::vector<Vertex*> vVertex;
1265 0 : CollectVertices(vVertex, grid, pElem);
1266 0 : vSubsetIndex.resize(vVertex.size());
1267 0 : for(size_t i = 0; i < vVertex.size(); ++i)
1268 0 : vSubsetIndex[i] = ish->get_subset_index(vVertex[i]);
1269 0 : }
1270 : if(dim == 2) {
1271 : std::vector<Edge*> vEdges;
1272 0 : CollectEdgesSorted(vEdges, grid, pElem);
1273 0 : vSubsetIndex.resize(vEdges.size());
1274 0 : for(size_t i = 0; i < vEdges.size(); ++i)
1275 0 : vSubsetIndex[i] = ish->get_subset_index(vEdges[i]);
1276 0 : }
1277 : if(dim == 3) {
1278 : std::vector<Face*> vFaces;
1279 0 : CollectFacesSorted(vFaces, grid, pElem);
1280 0 : vSubsetIndex.resize(vFaces.size());
1281 0 : for(size_t i = 0; i < vFaces.size(); ++i)
1282 0 : vSubsetIndex[i] = ish->get_subset_index(vFaces[i]);
1283 0 : }
1284 :
1285 : // update reference mapping
1286 0 : m_rMapping.update(vCornerCoords);
1287 :
1288 : // loop requested subset
1289 : typename std::map<int, std::vector<BF> >::iterator it;
1290 0 : for (it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); ++it)
1291 : {
1292 : // get subset index
1293 0 : const int bndIndex = (*it).first;
1294 :
1295 : // get vector of BF for element
1296 0 : std::vector<BF>& vBF = (*it).second;
1297 :
1298 : // clear vector
1299 : vBF.clear();
1300 :
1301 : // current number of bf
1302 : size_t curr_bf = 0;
1303 :
1304 : // loop subelements
1305 0 : for(size_t se = 0; se < numSubElem; ++se)
1306 : {
1307 : // skip inner sub elements
1308 0 : if(!m_vSubElem[se].isBndElem) continue;
1309 :
1310 : // loop sides of element
1311 0 : for(size_t side = 0; side < m_vSubElem[se].vElemBndSide.size(); ++side)
1312 : {
1313 : // get whole element bnd side
1314 0 : const int elemBndSide = m_vSubElem[se].vElemBndSide[side];
1315 :
1316 : // skip non boundary sides
1317 0 : if(elemBndSide == -1 || vSubsetIndex[elemBndSide] != bndIndex) continue;
1318 :
1319 : // number of corners of side
1320 0 : const int coOfSide = m_rRefElem.num(dim-1, elemBndSide, 0);
1321 :
1322 : // resize vector
1323 0 : vBF.resize(curr_bf + coOfSide);
1324 :
1325 : // loop corners
1326 0 : for(int co = 0; co < coOfSide; ++co)
1327 : {
1328 : // get current bf
1329 : BF& bf = vBF[curr_bf];
1330 :
1331 : // set node id == scv this bf belongs to
1332 0 : const int refNodeId = m_rRefElem.id(dim-1, elemBndSide, 0, co);
1333 0 : bf.nodeId = m_vSubElem[se].vDoFID[refNodeId];
1334 :
1335 : // Compute MidID for BF
1336 0 : ComputeBFMidID(m_rRefElem, elemBndSide, bf.vMidID, co);
1337 :
1338 : // copy corners of bf
1339 : CopyCornerByMidID<dim, maxMid>
1340 0 : (bf.vLocPos, bf.vMidID, m_vSubElem[se].vvLocMid, BF::numCo);
1341 : // CopyCornerByMidID<worldDim, maxMid>
1342 : // (bf.vGloPos, bf.vMidID, m_vSubElem[se].vvGloMid, BF::numCo);
1343 : // compute global corners of bf
1344 0 : for(size_t i = 0; i < bf.num_corners(); ++i)
1345 0 : m_rMapping.local_to_global(bf.vGloPos[i], bf.vLocPos[i]);
1346 :
1347 : // normal on scvf
1348 0 : traits::NormalOnSCVF(bf.Normal, bf.vGloPos, m_vSubElem[se].vvGloMid[0]);
1349 :
1350 : // compute local integration points
1351 0 : bf.vWeight = m_rSCVFQuadRule.weights();
1352 0 : ReferenceMapping<scvf_type, dim> map(bf.vLocPos);
1353 0 : for(size_t ip = 0; ip < m_rSCVFQuadRule.size(); ++ip)
1354 0 : map.local_to_global(bf.vLocalIP[ip], m_rSCVFQuadRule.point(ip));
1355 :
1356 : // compute global integration points
1357 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
1358 0 : m_rMapping.local_to_global(bf.vGlobalIP[ip], bf.vLocalIP[ip]);
1359 :
1360 : // compute volume
1361 0 : bf.Vol = VecTwoNorm(bf.Normal);
1362 :
1363 : // compute shapes and gradients
1364 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
1365 : {
1366 0 : m_rTrialSpace.shapes(&(bf.vvShape[ip][0]), bf.local_ip(ip));
1367 0 : m_rTrialSpace.grads(&(bf.vvLocalGrad[ip][0]), bf.local_ip(ip));
1368 :
1369 0 : m_rMapping.jacobian_transposed_inverse(bf.vJtInv[ip], bf.local_ip(ip));
1370 0 : bf.vDetJ[ip] = m_rMapping.sqrt_gram_det(bf.local_ip(ip));
1371 : }
1372 :
1373 : // compute global gradient
1374 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
1375 0 : for(size_t sh = 0 ; sh < bf.num_sh(); ++sh)
1376 : MatVecMult(bf.vvGlobalGrad[ip][sh],
1377 0 : bf.vJtInv[ip], bf.vvLocalGrad[ip][sh]);
1378 :
1379 : // increase curr_bf
1380 0 : ++curr_bf;
1381 : }
1382 : }
1383 : }
1384 : }
1385 0 : }
1386 :
1387 :
1388 : ////////////////////////////////////////////////////////////////////////////////
1389 : ////////////////////////////////////////////////////////////////////////////////
1390 : // FV Geometry (all order, FVHO) DIM FV
1391 : ////////////////////////////////////////////////////////////////////////////////
1392 : ////////////////////////////////////////////////////////////////////////////////
1393 :
1394 : template <int TWorldDim, int TDim>
1395 0 : DimFVGeometry<TWorldDim, TDim>::
1396 0 : DimFVGeometry() : m_pElem(NULL) {}
1397 :
1398 :
1399 :
1400 : template <int TWorldDim, int TDim>
1401 0 : void DimFVGeometry<TWorldDim, TDim>::
1402 : update_local(ReferenceObjectID roid, const LFEID& lfeID, size_t orderQuad)
1403 : {
1404 : // save setting we prepare the local data for
1405 0 : m_roid = roid;
1406 0 : m_lfeID = lfeID;
1407 0 : m_orderShape = m_lfeID.order();
1408 0 : m_quadOrderSCVF = orderQuad;
1409 0 : m_quadOrderSCV = orderQuad;
1410 :
1411 : // resize sub elements
1412 0 : m_numSubElem = (size_t) std::pow((double) m_orderShape, dim);
1413 0 : m_vSubElem.resize(m_numSubElem);
1414 :
1415 : // get the multi indices for the sub elements and the boundary flags
1416 0 : bool* vIsBndElem = new bool[m_numSubElem];
1417 0 : std::vector<std::vector<int> > vElemBndSide(m_numSubElem);
1418 0 : std::vector<std::vector<MathVector<dim,int> > > vMultiIndex(m_numSubElem);
1419 0 : std::vector<std::vector<size_t> > vIndex(m_numSubElem);
1420 0 : ComputeMultiIndicesOfSubElement<dim>(&vMultiIndex[0], &vIsBndElem[0],
1421 : &vElemBndSide[0], &vIndex[0], m_roid, m_orderShape);
1422 :
1423 : // get reference element
1424 : try{
1425 : const DimReferenceElement<dim>& rRefElem
1426 0 : = ReferenceElementProvider::get<dim>(m_roid);
1427 :
1428 : // directions of counting
1429 0 : MathVector<dim> direction[dim];
1430 0 : for(int i = 0; i < dim; ++i){direction[i] = 0.0; direction[i][i] = 1.0;}
1431 :
1432 0 : for(size_t se = 0; se < m_numSubElem; ++se)
1433 : {
1434 0 : for(size_t co = 0; co < vMultiIndex[se].size(); ++co)
1435 : {
1436 : // compute corners of sub elem in local coordinates
1437 : MathVector<dim> pos; pos = 0.0;
1438 0 : for(int i = 0; i < dim; ++i)
1439 : {
1440 0 : const number frac = vMultiIndex[se][co][i] / ((number)m_orderShape);
1441 : VecScaleAppend(pos, frac, direction[i]);
1442 : }
1443 : m_vSubElem[se].vvLocMid[0][co] = pos;
1444 :
1445 : // get multi index for corner
1446 0 : m_vSubElem[se].vDoFID[co] = vIndex[se][co];
1447 : }
1448 :
1449 : // remember if boundary element
1450 0 : m_vSubElem[se].isBndElem = vIsBndElem[se];
1451 :
1452 : // remember boundary sides
1453 0 : m_vSubElem[se].vElemBndSide = vElemBndSide[se];
1454 : }
1455 :
1456 0 : delete[] vIsBndElem;
1457 :
1458 : // compute mid points for all sub elements
1459 0 : for(size_t se = 0; se < m_numSubElem; ++se)
1460 : ComputeMidPoints<dim, DimReferenceElement<dim>, maxMid>
1461 0 : (rRefElem, m_vSubElem[se].vvLocMid[0], m_vSubElem[se].vvLocMid);
1462 :
1463 : // number of scvf/scv per subelem
1464 0 : m_numSCVFPerSubElem = rRefElem.num(1);
1465 0 : m_numSCVPerSubElem = rRefElem.num(0);
1466 :
1467 0 : m_numSCVF = m_numSCVFPerSubElem * m_numSubElem;
1468 0 : m_numSCV = m_numSCVPerSubElem * m_numSubElem;
1469 :
1470 0 : m_vSCVF.resize(m_numSCVF);
1471 0 : m_vSCV.resize(m_numSCV);
1472 :
1473 :
1474 : // get trial space
1475 : const LocalShapeFunctionSet<dim>& rTrialSpace =
1476 0 : LocalFiniteElementProvider::get<dim>(m_roid, m_lfeID);
1477 :
1478 : // request for quadrature rule
1479 : const ReferenceObjectID scvfRoid = scvf_type::REFERENCE_OBJECT_ID;
1480 : const QuadratureRule<dim-1>& rSCVFQuadRule
1481 0 : = QuadratureRuleProvider<dim-1>::get(scvfRoid, m_quadOrderSCVF);
1482 :
1483 0 : const int nipSCVF = rSCVFQuadRule.size();
1484 0 : m_numSCVFIP = m_numSCVF * nipSCVF;
1485 :
1486 : // set up local informations for SubControlVolumeFaces (scvf)
1487 : // each scvf is associated to one edge of the sub-element
1488 0 : for(size_t i = 0; i < num_scvf(); ++i)
1489 : {
1490 : // get corresponding subelement
1491 0 : const size_t se = i / m_numSCVFPerSubElem;
1492 0 : const size_t locSCVF = i % m_numSCVFPerSubElem;
1493 :
1494 : // this scvf separates the given nodes
1495 0 : const size_t locFrom = rRefElem.id(1, locSCVF, 0, 0);
1496 0 : const size_t locTo = rRefElem.id(1, locSCVF, 0, 1);
1497 :
1498 0 : m_vSCVF[i].From = m_vSubElem[se].vDoFID[locFrom];
1499 0 : m_vSCVF[i].To = m_vSubElem[se].vDoFID[locTo];
1500 :
1501 : // compute mid ids of the scvf
1502 0 : ComputeSCVFMidID(rRefElem, m_vSCVF[i].vMidID, locSCVF);
1503 :
1504 : // copy local corners of scvf
1505 : CopyCornerByMidID<dim, maxMid>
1506 0 : (m_vSCVF[i].vLocPos, m_vSCVF[i].vMidID, m_vSubElem[se].vvLocMid, SCVF::numCo);
1507 :
1508 : // compute integration points
1509 0 : m_vSCVF[i].vWeight = rSCVFQuadRule.weights();
1510 0 : m_vSCVF[i].nip = nipSCVF;
1511 :
1512 0 : m_vSCVF[i].vLocalIP.resize(nipSCVF);
1513 0 : m_vSCVF[i].vGlobalIP.resize(nipSCVF);
1514 :
1515 0 : m_vSCVF[i].vvShape.resize(nipSCVF);
1516 0 : m_vSCVF[i].vvLocalGrad.resize(nipSCVF);
1517 0 : m_vSCVF[i].vvGlobalGrad.resize(nipSCVF);
1518 0 : m_vSCVF[i].vJtInv.resize(nipSCVF);
1519 0 : m_vSCVF[i].vDetJ.resize(nipSCVF);
1520 0 : m_vSCVF[i].vDetJMap.resize(nipSCVF);
1521 :
1522 0 : m_vSCVF[i].nsh = rTrialSpace.num_sh();
1523 :
1524 0 : ReferenceMapping<scvf_type, dim> map(m_vSCVF[i].vLocPos);
1525 0 : for(size_t ip = 0; ip < rSCVFQuadRule.size(); ++ip)
1526 0 : map.local_to_global(m_vSCVF[i].vLocalIP[ip], rSCVFQuadRule.point(ip));
1527 : }
1528 :
1529 :
1530 : // request for quadrature rule
1531 : static const ReferenceObjectID scvRoid = scv_type::REFERENCE_OBJECT_ID;
1532 : const QuadratureRule<dim>& rSCVQuadRule
1533 0 : = QuadratureRuleProvider<dim>::get(scvRoid, m_quadOrderSCV);
1534 :
1535 0 : const int nipSCV = rSCVQuadRule.size();
1536 0 : m_numSCVIP = m_numSCV * nipSCV;
1537 :
1538 : // set up local informations for SubControlVolumes (scv)
1539 : // each scv is associated to one corner of the sub-element
1540 0 : for(size_t i = 0; i < num_scv(); ++i)
1541 : {
1542 : // get corresponding subelement
1543 0 : const size_t se = i / m_numSCVPerSubElem;
1544 0 : const size_t locSCV = i % m_numSCVPerSubElem;
1545 :
1546 : // store associated node
1547 0 : m_vSCV[i].nodeId = m_vSubElem[se].vDoFID[locSCV];;
1548 :
1549 : // compute mid ids scv
1550 0 : ComputeSCVMidID(rRefElem, m_vSCV[i].vMidID, locSCV);
1551 :
1552 : // copy local corners of scv
1553 : CopyCornerByMidID<dim, maxMid>
1554 0 : (m_vSCV[i].vLocPos, m_vSCV[i].vMidID, m_vSubElem[se].vvLocMid, m_vSCV[i].num_corners());
1555 :
1556 : // compute integration points
1557 0 : m_vSCV[i].vWeight = rSCVQuadRule.weights();
1558 0 : m_vSCV[i].nip = nipSCV;
1559 :
1560 0 : m_vSCV[i].vLocalIP.resize(nipSCV);
1561 0 : m_vSCV[i].vGlobalIP.resize(nipSCV);
1562 :
1563 0 : m_vSCV[i].vvShape.resize(nipSCV);
1564 0 : m_vSCV[i].vvLocalGrad.resize(nipSCV);
1565 0 : m_vSCV[i].vvGlobalGrad.resize(nipSCV);
1566 0 : m_vSCV[i].vJtInv.resize(nipSCV);
1567 0 : m_vSCV[i].vDetJ.resize(nipSCV);
1568 0 : m_vSCV[i].vDetJMap.resize(nipSCV);
1569 :
1570 0 : m_vSCV[i].nsh = rTrialSpace.num_sh();
1571 :
1572 0 : ReferenceMapping<scv_type, dim> map(m_vSCV[i].vLocPos);
1573 0 : for(size_t ip = 0; ip < rSCVQuadRule.size(); ++ip)
1574 0 : map.local_to_global(m_vSCV[i].vLocalIP[ip], rSCVQuadRule.point(ip));
1575 : }
1576 :
1577 : /////////////////////////
1578 : // Shapes and Derivatives
1579 : /////////////////////////
1580 :
1581 0 : m_nsh = rTrialSpace.num_sh();
1582 :
1583 0 : for(size_t i = 0; i < num_scvf(); ++i)
1584 0 : for(size_t ip = 0; ip < m_vSCVF[i].num_ip(); ++ip)
1585 : {
1586 0 : m_vSCVF[i].vvShape[ip].resize(m_vSCVF[i].nsh);
1587 0 : m_vSCVF[i].vvLocalGrad[ip].resize(m_vSCVF[i].nsh);
1588 0 : m_vSCVF[i].vvGlobalGrad[ip].resize(m_vSCVF[i].nsh);
1589 :
1590 0 : rTrialSpace.shapes(&(m_vSCVF[i].vvShape[ip][0]), m_vSCVF[i].local_ip(ip));
1591 0 : rTrialSpace.grads(&(m_vSCVF[i].vvLocalGrad[ip][0]), m_vSCVF[i].local_ip(ip));
1592 : }
1593 :
1594 0 : for(size_t i = 0; i < num_scv(); ++i)
1595 0 : for(size_t ip = 0; ip < m_vSCV[i].num_ip(); ++ip)
1596 : {
1597 0 : m_vSCV[i].vvShape[ip].resize(m_vSCV[i].nsh);
1598 0 : m_vSCV[i].vvLocalGrad[ip].resize(m_vSCV[i].nsh);
1599 0 : m_vSCV[i].vvGlobalGrad[ip].resize(m_vSCV[i].nsh);
1600 :
1601 0 : rTrialSpace.shapes(&(m_vSCV[i].vvShape[ip][0]), m_vSCV[i].local_ip(ip));
1602 0 : rTrialSpace.grads(&(m_vSCV[i].vvLocalGrad[ip][0]), m_vSCV[i].local_ip(ip));
1603 : }
1604 :
1605 : }
1606 0 : UG_CATCH_THROW("DimFV1Geometry: update failed.");
1607 :
1608 : // copy ip positions in a list for Sub Control Volumes Faces (SCVF)
1609 0 : m_vLocSCVF_IP.resize(m_numSCVFIP);
1610 0 : m_vGlobSCVF_IP.resize(m_numSCVFIP);
1611 : size_t allIP = 0;
1612 0 : for(size_t i = 0; i < num_scvf(); ++i)
1613 0 : for(size_t ip = 0; ip < m_vSCVF[i].num_ip(); ++ip)
1614 0 : m_vLocSCVF_IP[allIP++] = scvf(i).local_ip(ip);
1615 :
1616 : // copy ip positions in a list for Sub Control Volumes (SCV)
1617 0 : m_vLocSCV_IP.resize(m_numSCVIP);
1618 0 : m_vGlobSCV_IP.resize(m_numSCVIP);
1619 : allIP = 0;
1620 0 : for(size_t i = 0; i < num_scv(); ++i)
1621 0 : for(size_t ip = 0; ip < m_vSCV[i].num_ip(); ++ip)
1622 0 : m_vLocSCV_IP[allIP++] = scv(i).local_ip(ip);
1623 0 : }
1624 :
1625 :
1626 : /// update data for given element
1627 : template <int TWorldDim, int TDim>
1628 0 : void DimFVGeometry<TWorldDim, TDim>::
1629 : update(GridObject* pElem, const MathVector<worldDim>* vCornerCoords,
1630 : const LFEID& lfeID, size_t quadOrder,
1631 : const ISubsetHandler* ish)
1632 : {
1633 : // If already update for this element, do nothing
1634 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
1635 :
1636 : // get reference element type
1637 0 : ReferenceObjectID roid = (ReferenceObjectID)pElem->reference_object_id();
1638 :
1639 : // if already prepared for this roid, skip update of local values
1640 0 : if(m_roid != roid || lfeID != m_lfeID ||
1641 0 : (int)quadOrder != m_quadOrderSCVF || (int)quadOrder != m_quadOrderSCV)
1642 0 : update_local(roid, lfeID, quadOrder);
1643 :
1644 : // get reference element mapping
1645 : try{
1646 : DimReferenceMapping<dim, worldDim>& rMapping
1647 0 : = ReferenceMappingProvider::get<dim, worldDim>(roid);
1648 :
1649 : // update reference mapping
1650 0 : rMapping.update(vCornerCoords);
1651 :
1652 : // compute global informations for scvf
1653 0 : for(size_t i = 0; i < num_scvf(); ++i)
1654 : {
1655 : // map local corners of scvf to global
1656 0 : for(size_t co = 0; co < m_vSCVF[i].num_corners(); ++co)
1657 0 : rMapping.local_to_global(m_vSCVF[i].vGloPos[co], m_vSCVF[i].vLocPos[co]);
1658 :
1659 : // normal on scvf
1660 0 : traits::NormalOnSCVF(m_vSCVF[i].Normal, m_vSCVF[i].vGloPos, vCornerCoords);
1661 0 : VecNormalize(m_vSCVF[i].Normal, m_vSCVF[i].Normal);
1662 :
1663 : static const ReferenceObjectID scvfRoid = scvf_type::REFERENCE_OBJECT_ID;
1664 : const QuadratureRule<dim-1>& rSCVFQuadRule
1665 0 : = QuadratureRuleProvider<dim-1>::get(scvfRoid, m_quadOrderSCVF);
1666 0 : ReferenceMapping<scvf_type, worldDim> map(m_vSCVF[i].vGloPos);
1667 :
1668 0 : for(size_t ip = 0; ip < rSCVFQuadRule.size(); ++ip){
1669 0 : m_vSCVF[i].vDetJMap[ip] = map.sqrt_gram_det(rSCVFQuadRule.point(ip));
1670 0 : if(dim == 1) m_vSCVF[i].vDetJMap[ip] = 1.0;
1671 0 : map.local_to_global(m_vSCVF[i].vGlobalIP[ip], rSCVFQuadRule.point(ip));
1672 : }
1673 :
1674 0 : rMapping.jacobian_transposed_inverse(&m_vSCVF[i].vJtInv[0], &m_vSCVF[i].vLocalIP[0], m_vSCVF[i].num_ip());
1675 0 : rMapping.sqrt_gram_det(&m_vSCVF[i].vDetJ[0], &m_vSCVF[i].vLocalIP[0], m_vSCVF[i].num_ip());
1676 : }
1677 :
1678 : // compute size of scv
1679 0 : for(size_t i = 0; i < num_scv(); ++i)
1680 : {
1681 : // map local corners of scvf to global
1682 0 : rMapping.local_to_global(&m_vSCV[i].vGloPos[0], &m_vSCV[i].vLocPos[0], m_vSCV[i].num_corners());
1683 :
1684 0 : rMapping.jacobian_transposed_inverse(&m_vSCV[i].vJtInv[0], &m_vSCV[i].vLocalIP[0], m_vSCV[i].num_ip());
1685 0 : rMapping.sqrt_gram_det(&m_vSCV[i].vDetJ[0], &m_vSCV[i].vLocalIP[0], m_vSCV[i].num_ip());
1686 :
1687 :
1688 : static const ReferenceObjectID scvRoid = scv_type::REFERENCE_OBJECT_ID;
1689 : const QuadratureRule<dim>& rSCVQuadRule
1690 0 : = QuadratureRuleProvider<dim>::get(scvRoid, m_quadOrderSCV);
1691 :
1692 0 : ReferenceMapping<scv_type, worldDim> map(m_vSCV[i].vGloPos);
1693 0 : for(size_t ip = 0; ip < rSCVQuadRule.size(); ++ip){
1694 0 : m_vSCV[i].vDetJMap[ip] = map.sqrt_gram_det(rSCVQuadRule.point(ip));
1695 0 : map.local_to_global(m_vSCV[i].vGlobalIP[ip], rSCVQuadRule.point(ip));
1696 : }
1697 : }
1698 :
1699 : }
1700 0 : UG_CATCH_THROW("DimFVGeometry: update failed.");
1701 :
1702 : // compute global gradients
1703 0 : for(size_t i = 0; i < num_scvf(); ++i)
1704 0 : for(size_t ip = 0; ip < scvf(i).num_ip(); ++ip)
1705 0 : for(size_t sh = 0 ; sh < m_vSCVF[i].nsh; ++sh)
1706 : MatVecMult(m_vSCVF[i].vvGlobalGrad[ip][sh], m_vSCVF[i].vJtInv[ip], m_vSCVF[i].vvLocalGrad[ip][sh]);
1707 :
1708 0 : for(size_t i = 0; i < num_scv(); ++i)
1709 0 : for(size_t ip = 0; ip < scv(i).num_ip(); ++ip)
1710 0 : for(size_t sh = 0 ; sh < m_vSCV[i].nsh; ++sh)
1711 : MatVecMult(m_vSCV[i].vvGlobalGrad[ip][sh], m_vSCV[i].vJtInv[ip], m_vSCV[i].vvLocalGrad[ip][sh]);
1712 :
1713 : // Copy ip pos in list for SCVF
1714 : size_t allIP = 0;
1715 0 : for(size_t i = 0; i < num_scvf(); ++i)
1716 0 : for(size_t ip = 0; ip < scvf(i).num_ip(); ++ip)
1717 0 : m_vGlobSCVF_IP[allIP++] = scvf(i).global_ip(ip);
1718 :
1719 : allIP = 0;
1720 0 : for(size_t i = 0; i < num_scv(); ++i)
1721 0 : for(size_t ip = 0; ip < scv(i).num_ip(); ++ip)
1722 0 : m_vGlobSCV_IP[allIP++] = scv(i).global_ip(ip);
1723 :
1724 : // if no boundary subsets required, return
1725 0 : if(num_boundary_subsets() == 0 || ish == NULL) return;
1726 0 : else update_boundary_faces(pElem, vCornerCoords, ish);
1727 : }
1728 :
1729 : template <int TWorldDim, int TDim>
1730 0 : void DimFVGeometry<TWorldDim, TDim>::
1731 : update_boundary_faces(GridObject* pElem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
1732 : {
1733 : // get grid
1734 0 : Grid& grid = *(ish->grid());
1735 :
1736 : // vector of subset indices of side
1737 : std::vector<int> vSubsetIndex;
1738 :
1739 : // get subset indices for sides (i.e. edge in 2d, faces in 3d)
1740 : if(dim == 1) {
1741 : std::vector<Vertex*> vVertex;
1742 0 : CollectVertices(vVertex, grid, pElem);
1743 0 : vSubsetIndex.resize(vVertex.size());
1744 0 : for(size_t i = 0; i < vVertex.size(); ++i)
1745 0 : vSubsetIndex[i] = ish->get_subset_index(vVertex[i]);
1746 0 : }
1747 : if(dim == 2) {
1748 : std::vector<Edge*> vEdges;
1749 0 : CollectEdgesSorted(vEdges, grid, pElem);
1750 0 : vSubsetIndex.resize(vEdges.size());
1751 0 : for(size_t i = 0; i < vEdges.size(); ++i)
1752 0 : vSubsetIndex[i] = ish->get_subset_index(vEdges[i]);
1753 0 : }
1754 : if(dim == 3) {
1755 : std::vector<Face*> vFaces;
1756 0 : CollectFacesSorted(vFaces, grid, pElem);
1757 0 : vSubsetIndex.resize(vFaces.size());
1758 0 : for(size_t i = 0; i < vFaces.size(); ++i)
1759 0 : vSubsetIndex[i] = ish->get_subset_index(vFaces[i]);
1760 0 : }
1761 :
1762 : // get reference element mapping
1763 : try{
1764 : DimReferenceMapping<dim, worldDim>& rMapping
1765 0 : = ReferenceMappingProvider::get<dim, worldDim>(m_roid);
1766 :
1767 : const DimReferenceElement<dim>& rRefElem
1768 0 : = ReferenceElementProvider::get<dim>(m_roid);
1769 :
1770 : const ReferenceObjectID scvfRoid = scvf_type::REFERENCE_OBJECT_ID;
1771 : const QuadratureRule<dim-1>& rSCVFQuadRule
1772 0 : = QuadratureRuleProvider<dim-1>::get(scvfRoid, m_quadOrderSCVF);
1773 :
1774 : const LocalShapeFunctionSet<dim>& rTrialSpace =
1775 0 : LocalFiniteElementProvider::get<dim>(m_roid, LFEID(LFEID::LAGRANGE, dim, m_orderShape));
1776 :
1777 : // update reference mapping
1778 0 : rMapping.update(vCornerCoords);
1779 :
1780 : // loop requested subset
1781 : typename std::map<int, std::vector<BF> >::iterator it;
1782 0 : for (it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); ++it)
1783 : {
1784 : // get subset index
1785 0 : const int bndIndex = (*it).first;
1786 :
1787 : // get vector of BF for element
1788 0 : std::vector<BF>& vBF = (*it).second;
1789 :
1790 : // clear vector
1791 : vBF.clear();
1792 :
1793 : // current number of bf
1794 : size_t curr_bf = 0;
1795 :
1796 : // loop subelements
1797 0 : for(size_t se = 0; se < m_numSubElem; ++se)
1798 : {
1799 : // skip inner sub elements
1800 0 : if(!m_vSubElem[se].isBndElem) continue;
1801 :
1802 : // loop sides of element
1803 0 : for(size_t side = 0; side < m_vSubElem[se].vElemBndSide.size(); ++side)
1804 : {
1805 : // get whole element bnd side
1806 0 : const int elemBndSide = m_vSubElem[se].vElemBndSide[side];
1807 :
1808 : // skip non boundary sides
1809 0 : if(elemBndSide == -1 || vSubsetIndex[elemBndSide] != bndIndex) continue;
1810 :
1811 : // number of corners of side
1812 0 : const int coOfSide = rRefElem.num(dim-1, elemBndSide, 0);
1813 :
1814 : // compute global mids for vertices on this subelem
1815 : // (needed to calculate normal - at least if dim < worldDim)
1816 0 : for (size_t m = 0; m < (size_t) maxMid; ++m)
1817 0 : rMapping.local_to_global(m_vSubElem[se].vvGloMid[0][m], m_vSubElem[se].vvLocMid[0][m]);
1818 :
1819 : // resize vector
1820 0 : vBF.resize(curr_bf + coOfSide);
1821 :
1822 : // loop corners
1823 0 : for(int co = 0; co < coOfSide; ++co)
1824 : {
1825 : // get current bf
1826 : BF& bf = vBF[curr_bf];
1827 :
1828 : // set node id == scv this bf belongs to
1829 0 : const int refNodeId = rRefElem.id(dim-1, elemBndSide, 0, co);
1830 0 : bf.nodeId = m_vSubElem[se].vDoFID[refNodeId];
1831 :
1832 : // Compute MidID for BF
1833 0 : ComputeBFMidID(rRefElem, elemBndSide, bf.vMidID, co);
1834 :
1835 : // copy local corners of bf
1836 : CopyCornerByMidID<dim, maxMid>
1837 0 : (bf.vLocPos, bf.vMidID, m_vSubElem[se].vvLocMid, BF::numCo);
1838 : // CopyCornerByMidID<worldDim, maxMid>
1839 : // (bf.vGloPos, bf.vMidID, m_vSubElem[se].vvGloMid, BF::numCo);
1840 :
1841 : // compute global corners of bf
1842 0 : for(size_t i = 0; i < bf.num_corners(); ++i)
1843 0 : rMapping.local_to_global(bf.vGloPos[i], bf.vLocPos[i]);
1844 :
1845 : // normal on bf
1846 0 : traits::NormalOnSCVF(bf.Normal, bf.vGloPos, m_vSubElem[se].vvGloMid[0]);
1847 :
1848 : // compute volume
1849 0 : bf.Vol = VecTwoNorm(bf.Normal);
1850 :
1851 : // compute local integration points
1852 0 : bf.vWeight = rSCVFQuadRule.weights();
1853 0 : bf.nip = rSCVFQuadRule.size();
1854 0 : bf.vLocalIP.resize(bf.nip);
1855 0 : bf.vGlobalIP.resize(bf.nip);
1856 :
1857 0 : ReferenceMapping<scvf_type, dim> map(bf.vLocPos);
1858 0 : for(size_t ip = 0; ip < rSCVFQuadRule.size(); ++ip)
1859 0 : map.local_to_global(bf.vLocalIP[ip], rSCVFQuadRule.point(ip));
1860 :
1861 : // compute global integration points
1862 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
1863 0 : rMapping.local_to_global(bf.vGlobalIP[ip], bf.vLocalIP[ip]);
1864 :
1865 0 : bf.nsh = rTrialSpace.num_sh();
1866 0 : bf.vvShape.resize(bf.nip);
1867 0 : bf.vvLocalGrad.resize(bf.nip);
1868 0 : bf.vvGlobalGrad.resize(bf.nip);
1869 0 : bf.vJtInv.resize(bf.nip);
1870 0 : bf.vDetJ.resize(bf.nip);
1871 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
1872 : {
1873 0 : bf.vvShape[ip].resize(bf.nsh);
1874 0 : bf.vvLocalGrad[ip].resize(bf.nsh);
1875 0 : bf.vvGlobalGrad[ip].resize(bf.nsh);
1876 : }
1877 :
1878 : // compute shapes and gradients
1879 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
1880 : {
1881 0 : rTrialSpace.shapes(&(bf.vvShape[ip][0]), bf.local_ip(ip));
1882 0 : rTrialSpace.grads(&(bf.vvLocalGrad[ip][0]), bf.local_ip(ip));
1883 : }
1884 :
1885 0 : rMapping.jacobian_transposed_inverse(&bf.vJtInv[0], &bf.vLocalIP[0], bf.num_ip());
1886 0 : rMapping.sqrt_gram_det(&bf.vDetJ[0], &bf.vLocalIP[0], bf.num_ip());
1887 :
1888 : // compute global gradient
1889 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
1890 0 : for(size_t sh = 0 ; sh < bf.num_sh(); ++sh)
1891 : MatVecMult(bf.vvGlobalGrad[ip][sh],
1892 : bf.vJtInv[ip], bf.vvLocalGrad[ip][sh]);
1893 :
1894 : // increase curr_bf
1895 0 : ++curr_bf;
1896 : }
1897 : }
1898 : }
1899 : }
1900 :
1901 : }
1902 0 : UG_CATCH_THROW("DimFVGeometry: update failed.");
1903 0 : }
1904 :
1905 : //////////////////////
1906 : // FVGeometry
1907 : #ifdef UG_DIM_1
1908 : template class DimFVGeometry<1, 1>;
1909 : #endif
1910 :
1911 : #ifdef UG_DIM_2
1912 : template class FVGeometry<1, RegularEdge, 2>;
1913 : template class FVGeometry<1, Triangle, 2>;
1914 : template class FVGeometry<1, Quadrilateral, 2>;
1915 : template class FVGeometry<2, RegularEdge, 2>;
1916 : template class FVGeometry<2, Triangle, 2>;
1917 : template class FVGeometry<2, Quadrilateral, 2>;
1918 : template class FVGeometry<3, RegularEdge, 2>;
1919 : template class FVGeometry<3, Triangle, 2>;
1920 : template class FVGeometry<3, Quadrilateral, 2>;
1921 :
1922 : template class DimFVGeometry<2, 2>;
1923 : template class DimFVGeometry<2, 1>;
1924 : #endif
1925 :
1926 : #ifdef UG_DIM_3
1927 : template class FVGeometry<1, RegularEdge, 3>;
1928 : template class FVGeometry<1, Tetrahedron, 3>;
1929 : template class FVGeometry<1, Prism, 3>;
1930 : template class FVGeometry<1, Hexahedron, 3>;
1931 : template class FVGeometry<2, RegularEdge, 3>;
1932 : template class FVGeometry<2, Tetrahedron, 3>;
1933 : template class FVGeometry<2, Prism, 3>;
1934 : template class FVGeometry<2, Hexahedron, 3>;
1935 : template class FVGeometry<3, RegularEdge, 3>;
1936 : template class FVGeometry<3, Tetrahedron, 3>;
1937 : template class FVGeometry<3, Prism, 3>;
1938 : template class FVGeometry<3, Hexahedron, 3>;
1939 :
1940 : template class DimFVGeometry<3, 3>;
1941 : template class DimFVGeometry<3, 2>;
1942 : template class DimFVGeometry<3, 1>;
1943 : #endif
1944 :
1945 : } // end namespace ug
|
