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