Line data Source code
1 : /*
2 : * Copyright (c) 2010-2021: G-CSC, Goethe University Frankfurt
3 : * Authors: Andreas Vogel, Dmitrij Logashenko, Martin Stepniewski
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 "fv1_geom.h"
36 : #include "lib_disc/reference_element/reference_element.h"
37 : #include "lib_disc/quadrature/quadrature.h"
38 : #include "lib_algebra/common/operations_vec.h"
39 :
40 : namespace ug{
41 :
42 : DebugID DID_FV1_GEOM("FV1_GEOM");
43 : DebugID DID_REFERENCE_MAPPING("REFERENCE_MAPPING");
44 : DebugID DID_REFERENCE_MAPPING_GLOB_TO_LOC("REFERENCE_MAPPING_GLOB_TO_LOC");
45 : DebugID DID_ELEM_DISC_ASSEMBLE_UTIL("ELEM_DISC_ASSEMBLE_UTIL");
46 :
47 : /**
48 : * \tparam dim dimension of coordinates
49 : * \tparam TRefElem Reference element type
50 : * \tparam maxMid Maximum number of elements for all dimensions
51 : */
52 : template <int dim, typename TRefElem, int maxMid>
53 0 : static void ComputeMidPoints(const TRefElem& rRefElem,
54 : const MathVector<dim> vCorner[],
55 : MathVector<dim> vvMid[][maxMid])
56 : {
57 : // compute local midpoints for all geometric objects with 0 < d <= dim
58 0 : for(int d = 1; d <= dim; ++d)
59 : {
60 : // loop geometric objects of dimension d
61 0 : for(size_t i = 0; i < rRefElem.num(d); ++i)
62 : {
63 : // set first node
64 0 : const size_t coID0 = rRefElem.id(d, i, 0, 0);
65 0 : vvMid[d][i] = vCorner[coID0];
66 :
67 : // add corner coordinates of the corners of the geometric object
68 0 : for(size_t j = 1; j < rRefElem.num(d, i, 0); ++j)
69 : {
70 0 : const size_t coID = rRefElem.id(d, i, 0, j);
71 0 : vvMid[d][i] += vCorner[coID];
72 : }
73 :
74 : // scale for correct averaging
75 0 : vvMid[d][i] *= 1./(rRefElem.num(d, i, 0));
76 : }
77 : }
78 :
79 : /**
80 : * The octahedral reference element contains an implicit interior
81 : * substructure that is constructed by several geometric objects,
82 : * i.e. imaginary subfaces (8 triangles), subedges (2) and subvolumes (4 tetrahedra)
83 : * resulting from the division into 4 subtetrahedra alongside inner edge 3->1.
84 : * The dual fv1 geometry consists of the original hexahedral SCVs in each of the
85 : * 4 subtetrahedra. Therefore, add the following additional corresponding midpoints:
86 : */
87 0 : if(rRefElem.roid() == ROID_OCTAHEDRON)
88 : {
89 : typedef fv1_traits_ReferenceOctahedron traits;
90 :
91 0 : for(int d = 1; d <= dim; ++d)
92 : {
93 : // loop geometric objects of dimension d of the implicit interior substructure
94 0 : for(size_t i = 0; i < traits::substruct_num(d); ++i)
95 : {
96 : // set first node
97 0 : const size_t coID0 = traits::substruct_coID(d, i, 0);
98 0 : vvMid[d][rRefElem.num(d)+i] = vCorner[coID0];
99 :
100 : // add corner coordinates of the corners of the geometric object of the implicit interior substructure
101 0 : for(size_t j = 1; j < traits::substruct_num(d, i, 0); ++j)
102 : {
103 0 : const size_t coID = traits::substruct_coID(d, i, j);
104 0 : vvMid[d][rRefElem.num(d)+i] += vCorner[coID];
105 : }
106 :
107 : // scale for correct averaging
108 0 : vvMid[d][rRefElem.num(d)+i] *= 1./(traits::substruct_num(d, i, 0));
109 : }
110 : }
111 : }
112 0 : }
113 :
114 : /**
115 : * \param[in] i indicates that scvf corresponds to i'th edge of ref elem
116 : */
117 : template <typename TRefElem>
118 0 : static void ComputeSCVFMidID(const TRefElem& rRefElem,
119 : MidID vMidID[], int i)
120 : {
121 : static const int dim = TRefElem::dim;
122 :
123 0 : if (rRefElem.roid() != ROID_PYRAMID && rRefElem.roid() != ROID_OCTAHEDRON)
124 : {
125 : // start at edge midpoint
126 0 : vMidID[0] = MidID(1,i);
127 :
128 : // loop up dimension
129 : if(dim == 2)
130 : {
131 0 : vMidID[1] = MidID(dim, 0); // center of element
132 : }
133 : else if (dim == 3)
134 : {
135 0 : vMidID[1] = MidID(2, rRefElem.id(1, i, 2, 0)); // side 0
136 0 : vMidID[2] = MidID(dim, 0); // center of element
137 0 : vMidID[3] = MidID(2, rRefElem.id(1, i, 2, 1)); // side 1
138 : }
139 : }
140 : // pyramid here
141 0 : else if(rRefElem.roid() == ROID_PYRAMID)
142 : {
143 0 : fv1_traits_ReferencePyramid::scvf_mid_id((const ReferencePyramid&) rRefElem, vMidID, i);
144 : }
145 : // octahedron here (analogue to scvf ordering in pyramids but in top/bottom pairs)
146 : else if(rRefElem.roid() == ROID_OCTAHEDRON)
147 : {
148 0 : fv1_traits_ReferenceOctahedron::scvf_mid_id((const ReferenceOctahedron&) rRefElem, vMidID, i);
149 : }
150 0 : }
151 :
152 : /**
153 : * \param[in] i indicates that scv corresponds to i'th corner of ref elem
154 : */
155 : template <typename TRefElem>
156 0 : static void ComputeSCVMidID(const TRefElem& rRefElem,
157 : MidID vMidID[], int i)
158 : {
159 : static const int dim = TRefElem::dim;
160 :
161 0 : if (rRefElem.roid() != ROID_PYRAMID && rRefElem.roid() != ROID_OCTAHEDRON)
162 : {
163 : if(dim == 1)
164 : {
165 0 : vMidID[0] = MidID(0, i); // set node as corner of scv
166 0 : vMidID[1] = MidID(dim, 0); // center of element
167 : }
168 : else if(dim == 2)
169 : {
170 0 : vMidID[0] = MidID(0, i); // set node as corner of scv
171 0 : vMidID[1] = MidID(1, rRefElem.id(0, i, 1, 0)); // edge 1
172 0 : vMidID[2] = MidID(dim, 0); // center of element
173 0 : vMidID[3] = MidID(1, rRefElem.id(0, i, 1, 1)); // edge 2
174 : }
175 : else if(dim == 3)
176 : {
177 0 : vMidID[0] = MidID(0, i); // set node as corner of scv
178 0 : vMidID[1] = MidID(1, rRefElem.id(0, i, 1, 1)); // edge 1
179 0 : vMidID[2] = MidID(2, rRefElem.id(0, i, 2, 0)); // face 0
180 0 : vMidID[3] = MidID(1, rRefElem.id(0, i, 1, 0)); // edge 0
181 0 : vMidID[4] = MidID(1, rRefElem.id(0, i, 1, 2)); // edge 2
182 0 : vMidID[5] = MidID(2, rRefElem.id(0, i, 2, 2)); // face 2
183 0 : vMidID[6] = MidID(dim, 0); // center of element
184 0 : vMidID[7] = MidID(2, rRefElem.id(0, i, 2, 1)); // face 1
185 : }
186 : else {UG_THROW("Dimension higher than 3 not implemented.");}
187 : }
188 : // pyramid here
189 0 : else if (rRefElem.roid() == ROID_PYRAMID)
190 : {
191 0 : fv1_traits_ReferencePyramid::scv_mid_id((const ReferencePyramid&) rRefElem, vMidID, i);
192 : }
193 : // octahedron here (analogue to scv ordering in pyramids but in top/bottom pairs)
194 : else if(rRefElem.roid() == ROID_OCTAHEDRON)
195 : {
196 0 : fv1_traits_ReferenceOctahedron::scv_mid_id((const ReferenceOctahedron&) rRefElem, vMidID, i);
197 : }
198 0 : }
199 :
200 : /**
201 : * \param[in] i indicates that bf corresponds to i'th corner of ref elem
202 : */
203 : template <typename TRefElem>
204 : static void ComputeBFMidID(const TRefElem& rRefElem, int side,
205 : MidID vMidID[], int co)
206 : {
207 : static const int dim = TRefElem::dim;
208 :
209 : // number of corners of side
210 0 : const int coOfSide = rRefElem.num(dim-1, side, 0);
211 :
212 : // set mid ids
213 : if (dim == 1)
214 : {
215 0 : vMidID[0] = MidID(0, rRefElem.id(0, side, 0, co));
216 : }
217 : else if(dim == 2)
218 : {
219 0 : vMidID[co%2] = MidID(0, rRefElem.id(1, side, 0, co)); // corner of side
220 0 : vMidID[(co+1)%2] = MidID(1, side); // side midpoint
221 : }
222 : else if (dim == 3)
223 : {
224 0 : vMidID[0] = MidID(0, rRefElem.id(2, side, 0, co)); // corner of side
225 0 : vMidID[1] = MidID(1, rRefElem.id(2, side, 1, co)); // edge co
226 0 : vMidID[2] = MidID(2, side); // side midpoint
227 0 : vMidID[3] = MidID(1, rRefElem.id(2, side, 1, (co -1 + coOfSide)%coOfSide)); // edge co-1
228 : }
229 : }
230 :
231 : template <int dim, int maxMid>
232 0 : static void CopyCornerByMidID(MathVector<dim> vCorner[],
233 : const MidID vMidID[],
234 : MathVector<dim> vvMidPos[][maxMid],
235 : const size_t numCo)
236 : {
237 0 : for(size_t i = 0; i < numCo; ++i)
238 : {
239 0 : const size_t d = vMidID[i].dim;
240 0 : const size_t id = vMidID[i].id;
241 0 : vCorner[i] = vvMidPos[d][id];
242 : }
243 0 : }
244 :
245 : ////////////////////////////////////////////////////////////////////////////////
246 : // FV1 Geometry for Reference Element Type
247 : ////////////////////////////////////////////////////////////////////////////////
248 :
249 : template <typename TElem, int TWorldDim, bool TCondensed>
250 0 : FV1Geometry_gen<TElem, TWorldDim, TCondensed>::
251 : FV1Geometry_gen()
252 0 : : m_pElem(NULL), m_rRefElem(Provider<ref_elem_type>::get()),
253 0 : m_rTrialSpace(Provider<local_shape_fct_set_type>::get())
254 : {
255 0 : update_local_data();
256 0 : }
257 :
258 : template <typename TElem, int TWorldDim, bool TCondensed>
259 0 : void FV1Geometry_gen<TElem, TWorldDim, TCondensed>::
260 : update_local_data()
261 : {
262 : UG_DLOG(DID_FV1_GEOM, 2, ">>OCT_DISC_DEBUG: " << "fv1_geom.cpp: " << "update_local_data(): " << std::endl);
263 :
264 : // set corners of element as local centers of nodes
265 0 : for(size_t i = 0; i < m_rRefElem.num(0); ++i)
266 : m_vvLocMid[0][i] = m_rRefElem.corner(i);
267 :
268 : // compute local midpoints
269 0 : ComputeMidPoints<dim, ref_elem_type, maxMid>(m_rRefElem, m_vvLocMid[0], m_vvLocMid);
270 :
271 : // set up local information for SubControlVolumeFaces (scvf)
272 0 : for(size_t i = 0; i < num_scvf(); ++i)
273 : {
274 :
275 : // this scvf separates the given nodes
276 0 : m_vSCVF[i].From = traits::scvf_from_to(m_rRefElem, i, 0);
277 0 : m_vSCVF[i].To = traits::scvf_from_to(m_rRefElem, i, 1);
278 :
279 : // compute mid ids of the scvf
280 0 : ComputeSCVFMidID(m_rRefElem, m_vSCVF[i].vMidID, i);
281 :
282 : // copy local corners of scvf
283 0 : CopyCornerByMidID<dim, maxMid>(m_vSCVF[i].vLocPos, m_vSCVF[i].vMidID, m_vvLocMid, SCVF::numCo);
284 :
285 : // integration point
286 : if (! condensed_scvf_ips)
287 0 : AveragePositions(m_vSCVF[i].localIP, m_vSCVF[i].vLocPos, SCVF::numCo); // the center of the patch
288 : else
289 : m_vSCVF[i].localIP = m_vSCVF[i].vLocPos[0]; // the midpoint of the edge
290 : }
291 :
292 : // set up local informations for SubControlVolumes (scv)
293 : // each scv is associated to one corner of the element
294 0 : for(size_t i = 0; i < num_scv(); ++i)
295 : {
296 : // store associated node
297 0 : m_vSCV[i].nodeId = traits::scv_node_id (m_rRefElem, i); // for 'classical elements', m_vSCV[i].nodeId = i
298 :
299 : // compute mid ids scv
300 0 : ComputeSCVMidID(m_rRefElem, m_vSCV[i].midId, i);
301 :
302 : // copy local corners of scv
303 0 : CopyCornerByMidID<dim, maxMid>(m_vSCV[i].vLocPos, m_vSCV[i].midId, m_vvLocMid, m_vSCV[i].num_corners());
304 : }
305 :
306 : // compute Shapes and Derivatives
307 0 : for(size_t i = 0; i < num_scvf(); ++i)
308 : {
309 0 : m_rTrialSpace.shapes(&(m_vSCVF[i].vShape[0]), m_vSCVF[i].local_ip());
310 0 : m_rTrialSpace.grads(&(m_vSCVF[i].vLocalGrad[0]), m_vSCVF[i].local_ip());
311 : }
312 :
313 0 : for(size_t i = 0; i < num_scv(); ++i)
314 : {
315 0 : m_rTrialSpace.shapes(&(m_vSCV[i].vShape[0]), m_vSCV[i].local_ip());
316 0 : m_rTrialSpace.grads(&(m_vSCV[i].vLocalGrad[0]), m_vSCV[i].local_ip());
317 : }
318 :
319 : // copy ip positions in a list for Sub Control Volumes Faces (SCVF)
320 0 : for(size_t i = 0; i < num_scvf(); ++i)
321 : m_vLocSCVF_IP[i] = scvf(i).local_ip();
322 :
323 : if(ref_elem_type::REFERENCE_OBJECT_ID == ROID_PYRAMID || ref_elem_type::REFERENCE_OBJECT_ID == ROID_OCTAHEDRON)
324 0 : for(size_t i = 0; i < num_scv(); ++i)
325 : m_vLocSCV_IP[i] = scv(i).local_ip();
326 0 : }
327 :
328 : /// update data for given element
329 : template <typename TElem, int TWorldDim, bool TCondensed>
330 0 : void FV1Geometry_gen<TElem, TWorldDim, TCondensed>::
331 : update(GridObject* elem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
332 : {
333 : UG_ASSERT(dynamic_cast<TElem*>(elem) != NULL, "Wrong element type.");
334 : TElem* pElem = static_cast<TElem*>(elem);
335 :
336 : // if already update for this element, do nothing
337 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
338 :
339 : // remember global position of nodes
340 0 : for(size_t i = 0; i < m_rRefElem.num(0); ++i)
341 0 : m_vvGloMid[0][i] = vCornerCoords[i];
342 :
343 : // compute global midpoints
344 0 : ComputeMidPoints<worldDim, ref_elem_type, maxMid>(m_rRefElem, m_vvGloMid[0], m_vvGloMid);
345 :
346 : UG_DLOG(DID_FV1_GEOM, 2, ">>OCT_DISC_DEBUG: " << "fv1_geom.cpp: " << "update(): " << "ComputeMidPoints(): " << std::endl);
347 : for(size_t i = 0; i < m_rRefElem.num(1)+2; ++i)
348 : {
349 : UG_DLOG(DID_FV1_GEOM, 2, " Edge midpoint " << i << ": " << m_vvGloMid[1][i] << std::endl);
350 : }
351 : for(size_t i = 0; i < m_rRefElem.num(2)+8; ++i)
352 : {
353 : UG_DLOG(DID_FV1_GEOM, 2, " Face midpoint " << i << ": " << ((dim >= 2) ? m_vvGloMid[2][i] : 0) << std::endl);
354 : }
355 : for(size_t i = 0; i < m_rRefElem.num(3)+4; ++i)
356 : {
357 : UG_DLOG(DID_FV1_GEOM, 2, " Volume midpoint " << i << ": " << ((dim >= 3) ? m_vvGloMid[3][i] : 0) << std::endl);
358 : }
359 :
360 : // compute global informations for scvf
361 : UG_DLOG(DID_FV1_GEOM, 2, ">>OCT_DISC_DEBUG: " << "fv1_geom.cpp: " << "update(): " << "scvf global info: " << std::endl);
362 0 : for(size_t i = 0; i < num_scvf(); ++i)
363 : {
364 : // copy local corners of scvf
365 0 : CopyCornerByMidID<worldDim, maxMid>(m_vSCVF[i].vGloPos, m_vSCVF[i].vMidID, m_vvGloMid, SCVF::numCo);
366 :
367 : // integration point
368 : if (! condensed_scvf_ips)
369 0 : AveragePositions(m_vSCVF[i].globalIP, m_vSCVF[i].vGloPos, SCVF::numCo); // the center of the patch
370 : else
371 : m_vSCVF[i].globalIP = m_vSCVF[i].vGloPos[0]; // the midpoint of the edge
372 :
373 : // normal on scvf
374 0 : traits::NormalOnSCVF(m_vSCVF[i].Normal, m_vSCVF[i].vGloPos, m_vvGloMid[0]);
375 : UG_DLOG(DID_FV1_GEOM, 2, " scvf # " << i << ": " << "m_vSCVF[i].globalIP: " << m_vSCVF[i].globalIP << "; m_vSCVF[i].localIP: " << m_vSCVF[i].localIP << "; \t \t m_vSCVF[i].Normal: " << m_vSCVF[i].Normal << "; m_vSCVF[i].NormalSize: " << VecLength(m_vSCVF[i].Normal) << std::endl);
376 : }
377 :
378 : // compute size of scv
379 : UG_DLOG(DID_FV1_GEOM, 2, ">>OCT_DISC_DEBUG: " << "fv1_geom.cpp: " << "update(): " << "scv global info: " << std::endl);
380 0 : for(size_t i = 0; i < num_scv(); ++i)
381 : {
382 : // copy global corners
383 0 : CopyCornerByMidID<worldDim, maxMid>(m_vSCV[i].vGloPos, m_vSCV[i].midId, m_vvGloMid, m_vSCV[i].num_corners());
384 :
385 : // compute volume of scv
386 0 : m_vSCV[i].Vol = ElementSize<scv_type, worldDim>(m_vSCV[i].vGloPos);
387 :
388 : /*
389 : * Only for debug purposes testing octahedral FV1 discretization
390 : *
391 : MathVector<dim> baryCenter(0.0);
392 : for(size_t j = 0; j < 8; ++j)
393 : {
394 : baryCenter += m_vSCV[i].vLocPos[j];
395 : }
396 :
397 : baryCenter *= 1.0/8.0;
398 : *
399 : */
400 :
401 : UG_DLOG(DID_FV1_GEOM, 2, " scv # " << i << ": " << "m_vSCV[i].vGloPos: " << m_vSCV[i].vGloPos[0] << "; m_vSCV[i].vLocPos: " << m_vSCV[i].vLocPos[0] << "; m_vSCV[i].Vol: " << m_vSCV[i].Vol << /*"; baryCenter: " << baryCenter <<*/ std::endl);
402 : }
403 :
404 : // Shapes and Derivatives
405 0 : m_mapping.update(vCornerCoords);
406 :
407 : // if mapping is linear, compute jacobian only once and copy
408 : if(ReferenceMapping<ref_elem_type, worldDim>::isLinear)
409 : {
410 : MathMatrix<worldDim,dim> JtInv;
411 0 : m_mapping.jacobian_transposed_inverse(JtInv, m_vSCVF[0].local_ip());
412 0 : const number detJ = m_mapping.sqrt_gram_det(m_vSCVF[0].local_ip());
413 :
414 0 : for(size_t i = 0; i < num_scvf(); ++i)
415 : {
416 : m_vSCVF[i].JtInv = JtInv;
417 0 : m_vSCVF[i].detj = detJ;
418 : }
419 :
420 0 : for(size_t i = 0; i < num_scv(); ++i)
421 : {
422 : m_vSCV[i].JtInv = JtInv;
423 0 : m_vSCV[i].detj = detJ;
424 : }
425 : }
426 : // else compute jacobian for each integration point
427 : else
428 : {
429 : UG_DLOG(DID_FV1_GEOM, 2, ">>OCT_DISC_DEBUG: " << "fv1_geom.cpp: " << "update(): num_scvf: " << num_scvf() << "; num_scv(): " << num_scv() << std::endl);
430 0 : for(size_t i = 0; i < num_scvf(); ++i)
431 : {
432 : UG_DLOG(DID_FV1_GEOM, 2, ">>OCT_DISC_DEBUG: " << "fv1_geom.cpp: " << "update():" << "Jacobian for ip in SCVF # " << i << "; local_ip(): " << m_vSCVF[i].local_ip() << "; global_ip(): " << m_vSCVF[i].global_ip() << std::endl);
433 0 : m_mapping.jacobian_transposed_inverse(m_vSCVF[i].JtInv, m_vSCVF[i].local_ip());
434 0 : m_vSCVF[i].detj = m_mapping.sqrt_gram_det(m_vSCVF[i].local_ip());
435 : }
436 0 : for(size_t i = 0; i < num_scv(); ++i)
437 : {
438 : UG_DLOG(DID_FV1_GEOM, 2, ">>OCT_DISC_DEBUG: " << "fv1_geom.cpp: " << "update():" << "Jacobian for ip in SCV # " << i << "; local_ip(): " << m_vSCV[i].local_ip() << "; global_ip(): " << m_vSCV[i].global_ip() << std::endl);
439 0 : m_mapping.jacobian_transposed_inverse(m_vSCV[i].JtInv, m_vSCV[i].local_ip());
440 0 : m_vSCV[i].detj = m_mapping.sqrt_gram_det(m_vSCV[i].local_ip());
441 : }
442 : }
443 :
444 : // compute global gradients
445 0 : for(size_t i = 0; i < num_scvf(); ++i)
446 0 : for(size_t sh = 0 ; sh < scvf(i).num_sh(); ++sh)
447 0 : MatVecMult(m_vSCVF[i].vGlobalGrad[sh], m_vSCVF[i].JtInv, m_vSCVF[i].vLocalGrad[sh]);
448 :
449 0 : for(size_t i = 0; i < num_scv(); ++i)
450 0 : for(size_t sh = 0 ; sh < scv(i).num_sh(); ++sh)
451 0 : MatVecMult(m_vSCV[i].vGlobalGrad[sh], m_vSCV[i].JtInv, m_vSCV[i].vLocalGrad[sh]);
452 :
453 : // Copy ip pos in list for SCVF
454 0 : for(size_t i = 0; i < num_scvf(); ++i)
455 : m_vGlobSCVF_IP[i] = scvf(i).global_ip();
456 :
457 : if(ref_elem_type::REFERENCE_OBJECT_ID == ROID_PYRAMID || ref_elem_type::REFERENCE_OBJECT_ID == ROID_OCTAHEDRON)
458 0 : for(size_t i = 0; i < num_scv(); ++i)
459 : m_vGlobSCV_IP[i] = scv(i).global_ip();
460 :
461 : // if no boundary subsets required, return
462 0 : if(num_boundary_subsets() == 0 || ish == NULL) return;
463 0 : else update_boundary_faces(pElem, vCornerCoords, ish);
464 : }
465 :
466 : template <typename TElem, int TWorldDim, bool TCondensed>
467 0 : void FV1Geometry_gen<TElem, TWorldDim, TCondensed>::
468 : update_boundary_faces(GridObject* elem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
469 : {
470 : UG_ASSERT(dynamic_cast<TElem*>(elem) != NULL, "Wrong element type.");
471 : TElem* pElem = static_cast<TElem*>(elem);
472 :
473 : // get grid
474 0 : Grid& grid = *(ish->grid());
475 :
476 : // vector of subset indices of side
477 : std::vector<int> vSubsetIndex;
478 :
479 : // get subset indices for sides (i.e. edge in 2d, faces in 3d)
480 : if(dim == 1) {
481 : std::vector<Vertex*> vVertex;
482 0 : CollectVertices(vVertex, grid, pElem);
483 0 : vSubsetIndex.resize(vVertex.size());
484 0 : for(size_t i = 0; i < vVertex.size(); ++i)
485 0 : vSubsetIndex[i] = ish->get_subset_index(vVertex[i]);
486 0 : }
487 : if(dim == 2) {
488 : std::vector<Edge*> vEdges;
489 0 : CollectEdgesSorted(vEdges, grid, pElem);
490 0 : vSubsetIndex.resize(vEdges.size());
491 0 : for(size_t i = 0; i < vEdges.size(); ++i)
492 0 : vSubsetIndex[i] = ish->get_subset_index(vEdges[i]);
493 0 : }
494 : if(dim == 3) {
495 : std::vector<Face*> vFaces;
496 0 : CollectFacesSorted(vFaces, grid, pElem);
497 0 : vSubsetIndex.resize(vFaces.size());
498 0 : for(size_t i = 0; i < vFaces.size(); ++i)
499 0 : vSubsetIndex[i] = ish->get_subset_index(vFaces[i]);
500 0 : }
501 :
502 : // loop requested subset
503 : typename std::map<int, std::vector<BF> >::iterator it;
504 0 : for (it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); ++it)
505 : {
506 : // get subset index
507 0 : const int bndIndex = (*it).first;
508 :
509 : // get vector of BF for element
510 0 : std::vector<BF>& vBF = (*it).second;
511 :
512 : // clear vector
513 : vBF.clear();
514 :
515 : // current number of bf
516 : size_t curr_bf = 0;
517 :
518 : // loop sides of element
519 0 : for(size_t side = 0; side < vSubsetIndex.size(); ++side)
520 : {
521 : // skip non boundary sides
522 0 : if(vSubsetIndex[side] != bndIndex) continue;
523 :
524 : // number of corners of side
525 0 : const int coOfSide = m_rRefElem.num(dim-1, side, 0);
526 :
527 : // resize vector
528 0 : vBF.resize(curr_bf + coOfSide);
529 :
530 : // loop corners
531 0 : for(int co = 0; co < coOfSide; ++co)
532 : {
533 : // get current bf
534 : BF& bf = vBF[curr_bf];
535 :
536 : // set node id == scv this bf belongs to
537 0 : bf.nodeId = m_rRefElem.id(dim-1, side, 0, co);
538 :
539 : // Compute MidID for BF
540 0 : ComputeBFMidID(m_rRefElem, side, bf.vMidID, co);
541 :
542 : // copy corners of bf
543 0 : CopyCornerByMidID<dim, maxMid>(bf.vLocPos, bf.vMidID, m_vvLocMid, BF::numCo);
544 0 : CopyCornerByMidID<worldDim, maxMid>(bf.vGloPos, bf.vMidID, m_vvGloMid, BF::numCo);
545 :
546 : // integration point
547 0 : AveragePositions(bf.localIP, bf.vLocPos, BF::numCo);
548 0 : AveragePositions(bf.globalIP, bf.vGloPos, BF::numCo);
549 :
550 : // normal on scvf
551 0 : traits::NormalOnBF(bf.Normal, bf.vGloPos, m_vvGloMid[0]);
552 :
553 : // compute volume
554 0 : bf.Vol = VecTwoNorm(bf.Normal);
555 :
556 0 : m_rTrialSpace.shapes(&(bf.vShape[0]), bf.localIP);
557 0 : m_rTrialSpace.grads(&(bf.vLocalGrad[0]), bf.localIP);
558 :
559 0 : m_mapping.jacobian_transposed_inverse(bf.JtInv, bf.localIP);
560 0 : bf.detj = m_mapping.sqrt_gram_det(bf.localIP);
561 :
562 0 : for(size_t sh = 0 ; sh < bf.num_sh(); ++sh)
563 0 : MatVecMult(bf.vGlobalGrad[sh], bf.JtInv, bf.vLocalGrad[sh]);
564 :
565 : // increase curr_bf
566 0 : ++curr_bf;
567 : }
568 : }
569 : }
570 0 : }
571 :
572 : ////////////////////////////////////////////////////////////////////////////////
573 : // Dim-dependent Finite Volume Geometry
574 : ////////////////////////////////////////////////////////////////////////////////
575 :
576 : template <int TDim, int TWorldDim>
577 0 : void DimFV1Geometry<TDim, TWorldDim>::
578 : update_local(ReferenceObjectID roid)
579 : {
580 0 : m_roid = roid;
581 :
582 : // get reference element
583 : try
584 : {
585 : const DimReferenceElement<dim>& rRefElem
586 : = ReferenceElementProvider::get<dim>(m_roid);
587 :
588 : // set corners of element as local centers of nodes
589 0 : for(size_t i = 0; i < rRefElem.num(0); ++i)
590 : m_vvLocMid[0][i] = rRefElem.corner(i);
591 :
592 : // compute local midpoints
593 : ComputeMidPoints<dim, DimReferenceElement<dim>, maxMid>
594 0 : (rRefElem, m_vvLocMid[0], m_vvLocMid);
595 :
596 : // set number of scvf / scv of this roid
597 0 : traits::dim_get_num_SCV_and_SCVF(rRefElem, m_roid, m_numSCV, m_numSCVF);
598 :
599 : // set up local informations for SubControlVolumeFaces (scvf)
600 : // each scvf is associated to one edge of the element
601 0 : for(size_t i = 0; i < num_scvf(); ++i)
602 : {
603 : // this scvf separates the given nodes
604 0 : traits::get_dim_scvf_from_to(rRefElem, m_roid, i, m_vSCVF[i].From, m_vSCVF[i].To);
605 :
606 : // compute mid ids of the scvf
607 0 : ComputeSCVFMidID(rRefElem, m_vSCVF[i].vMidID, i);
608 :
609 : // copy local corners of scvf
610 0 : CopyCornerByMidID<dim, maxMid>(m_vSCVF[i].vLocPos, m_vSCVF[i].vMidID, m_vvLocMid, SCVF::numCo);
611 :
612 : // integration point
613 0 : AveragePositions(m_vSCVF[i].localIP, m_vSCVF[i].vLocPos, SCVF::numCo);
614 : }
615 :
616 : // set up local informations for SubControlVolumes (scv)
617 : // each scv is associated to one corner of the element
618 0 : for(size_t i = 0; i < num_scv(); ++i)
619 : {
620 : // store associated node
621 0 : m_vSCV[i].nodeId = traits::dim_scv_node_id(rRefElem, m_roid, i); // for 'classical elements', m_vSCV[i].nodeId = i
622 :
623 : // compute mid ids scv
624 0 : ComputeSCVMidID(rRefElem, m_vSCV[i].vMidID, i);
625 :
626 : // copy local corners of scv
627 0 : CopyCornerByMidID<dim, maxMid>(m_vSCV[i].vLocPos, m_vSCV[i].vMidID, m_vvLocMid, m_vSCV[i].num_corners());
628 : }
629 :
630 : /////////////////////////
631 : // Shapes and Derivatives
632 : /////////////////////////
633 :
634 : const LocalShapeFunctionSet<dim>& TrialSpace =
635 0 : LocalFiniteElementProvider::get<dim>(m_roid, LFEID(LFEID::LAGRANGE, dim, 1));
636 :
637 0 : m_nsh = TrialSpace.num_sh();
638 :
639 0 : for(size_t i = 0; i < num_scvf(); ++i)
640 : {
641 0 : m_vSCVF[i].numSH = TrialSpace.num_sh();
642 0 : TrialSpace.shapes(&(m_vSCVF[i].vShape[0]), m_vSCVF[i].localIP);
643 0 : TrialSpace.grads(&(m_vSCVF[i].vLocalGrad[0]), m_vSCVF[i].localIP);
644 : }
645 :
646 0 : for(size_t i = 0; i < num_scv(); ++i)
647 : {
648 0 : m_vSCV[i].numSH = TrialSpace.num_sh();
649 0 : TrialSpace.shapes(&(m_vSCV[i].vShape[0]), m_vSCV[i].vLocPos[0]);
650 0 : TrialSpace.grads(&(m_vSCV[i].vLocalGrad[0]), m_vSCV[i].vLocPos[0]);
651 : }
652 :
653 : }
654 0 : UG_CATCH_THROW("DimFV1Geometry: update failed.");
655 :
656 : // copy ip positions in a list for Sub Control Volumes Faces (SCVF)
657 0 : for(size_t i = 0; i < num_scvf(); ++i)
658 : m_vLocSCVF_IP[i] = scvf(i).local_ip();
659 :
660 0 : if(roid == ROID_PYRAMID || roid == ROID_OCTAHEDRON)
661 0 : for(size_t i = 0; i < num_scv(); ++i)
662 : m_vLocSCV_IP[i] = scv(i).local_ip();
663 0 : }
664 :
665 :
666 : /// update data for given element
667 : template <int TDim, int TWorldDim>
668 0 : void DimFV1Geometry<TDim, TWorldDim>::
669 : update(GridObject* pElem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
670 : {
671 : // If already update for this element, do nothing
672 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
673 :
674 : // refresh local data, if different roid given
675 0 : if(m_roid != pElem->reference_object_id()) // remember new roid and update local data
676 0 : update_local ((ReferenceObjectID) pElem->reference_object_id());
677 :
678 : // get reference element
679 : try{
680 : const DimReferenceElement<dim>& rRefElem
681 0 : = ReferenceElementProvider::get<dim>(m_roid);
682 :
683 : // remember global position of nodes
684 0 : for(size_t i = 0; i < rRefElem.num(0); ++i)
685 0 : m_vvGloMid[0][i] = vCornerCoords[i];
686 :
687 : // compute local midpoints
688 0 : ComputeMidPoints<worldDim, DimReferenceElement<dim>, maxMid>(rRefElem, m_vvGloMid[0], m_vvGloMid);
689 :
690 : // compute global informations for scvf
691 0 : for(size_t i = 0; i < num_scvf(); ++i)
692 : {
693 : // copy local corners of scvf
694 0 : CopyCornerByMidID<worldDim, maxMid>(m_vSCVF[i].vGloPos, m_vSCVF[i].vMidID, m_vvGloMid, SCVF::numCo);
695 :
696 : // integration point
697 0 : AveragePositions(m_vSCVF[i].globalIP, m_vSCVF[i].vGloPos, SCVF::numCo);
698 :
699 : // normal on scvf
700 0 : traits::NormalOnSCVF(m_vSCVF[i].Normal, m_vSCVF[i].vGloPos, m_vvGloMid[0]);
701 : }
702 :
703 : // compute size of scv
704 0 : for(size_t i = 0; i < num_scv(); ++i)
705 : {
706 : // copy global corners
707 0 : CopyCornerByMidID<worldDim, maxMid>(m_vSCV[i].vGloPos, m_vSCV[i].vMidID, m_vvGloMid, m_vSCV[i].num_corners());
708 :
709 : // compute volume of scv
710 0 : m_vSCV[i].Vol = ElementSize<scv_type, worldDim>(m_vSCV[i].vGloPos);
711 : }
712 :
713 : // get reference mapping
714 0 : DimReferenceMapping<dim, worldDim>& rMapping = ReferenceMappingProvider::get<dim, worldDim>(m_roid);
715 0 : rMapping.update(vCornerCoords);
716 :
717 : //\todo compute with on virt. call
718 : // compute jacobian for linear mapping
719 0 : if(rMapping.is_linear())
720 : {
721 : MathMatrix<worldDim,dim> JtInv;
722 0 : rMapping.jacobian_transposed_inverse(JtInv, m_vSCVF[0].local_ip());
723 0 : const number detJ = rMapping.sqrt_gram_det(m_vSCVF[0].local_ip());
724 :
725 0 : for(size_t i = 0; i < num_scvf(); ++i)
726 : {
727 : m_vSCVF[i].JtInv = JtInv;
728 0 : m_vSCVF[i].detj = detJ;
729 : }
730 :
731 0 : for(size_t i = 0; i < num_scv(); ++i)
732 : {
733 : m_vSCV[i].JtInv = JtInv;
734 0 : m_vSCV[i].detj = detJ;
735 : }
736 : }
737 : // else compute jacobian for each integration point
738 : else
739 : {
740 0 : for(size_t i = 0; i < num_scvf(); ++i)
741 : {
742 0 : rMapping.jacobian_transposed_inverse(m_vSCVF[i].JtInv, m_vSCVF[i].local_ip());
743 0 : m_vSCVF[i].detj = rMapping.sqrt_gram_det(m_vSCVF[i].local_ip());
744 : }
745 0 : for(size_t i = 0; i < num_scv(); ++i)
746 : {
747 0 : rMapping.jacobian_transposed_inverse(m_vSCV[i].JtInv, m_vSCV[i].local_ip());
748 0 : m_vSCV[i].detj = rMapping.sqrt_gram_det(m_vSCV[i].local_ip());
749 : }
750 : }
751 :
752 : // compute global gradients
753 0 : for(size_t i = 0; i < num_scvf(); ++i)
754 0 : for(size_t sh = 0; sh < scvf(i).num_sh(); ++sh)
755 0 : MatVecMult(m_vSCVF[i].vGlobalGrad[sh], m_vSCVF[i].JtInv, m_vSCVF[i].vLocalGrad[sh]);
756 :
757 0 : for(size_t i = 0; i < num_scv(); ++i)
758 0 : for(size_t sh = 0; sh < scv(i).num_sh(); ++sh)
759 0 : MatVecMult(m_vSCV[i].vGlobalGrad[sh], m_vSCV[i].JtInv, m_vSCV[i].vLocalGrad[sh]);
760 :
761 : // copy ip points in list (SCVF)
762 0 : for(size_t i = 0; i < num_scvf(); ++i)
763 : m_vGlobSCVF_IP[i] = scvf(i).global_ip();
764 :
765 0 : if(m_roid == ROID_PYRAMID || m_roid == ROID_OCTAHEDRON)
766 0 : for(size_t i = 0; i < num_scv(); ++i)
767 : m_vGlobSCV_IP[i] = scv(i).global_ip();
768 :
769 : }
770 0 : UG_CATCH_THROW("DimFV1Geometry: update failed.");
771 :
772 : // if no boundary subsets required, return
773 0 : if(num_boundary_subsets() == 0 || ish == NULL) return;
774 0 : else update_boundary_faces(pElem, vCornerCoords, ish);
775 : }
776 :
777 : template <int TDim, int TWorldDim>
778 0 : void DimFV1Geometry<TDim, TWorldDim>::
779 : update_boundary_faces(GridObject* pElem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
780 : {
781 : // get grid
782 0 : Grid& grid = *(ish->grid());
783 :
784 : // vector of subset indices of side
785 : std::vector<int> vSubsetIndex;
786 :
787 : // get subset indices for sides (i.e. edge in 2d, faces in 3d)
788 : if(dim == 1) {
789 : std::vector<Vertex*> vVertex;
790 0 : CollectVertices(vVertex, grid, pElem);
791 0 : vSubsetIndex.resize(vVertex.size());
792 0 : for(size_t i = 0; i < vVertex.size(); ++i)
793 0 : vSubsetIndex[i] = ish->get_subset_index(vVertex[i]);
794 0 : }
795 : if(dim == 2) {
796 : std::vector<Edge*> vEdges;
797 0 : CollectEdgesSorted(vEdges, grid, pElem);
798 0 : vSubsetIndex.resize(vEdges.size());
799 0 : for(size_t i = 0; i < vEdges.size(); ++i)
800 0 : vSubsetIndex[i] = ish->get_subset_index(vEdges[i]);
801 0 : }
802 : if(dim == 3) {
803 : std::vector<Face*> vFaces;
804 0 : CollectFacesSorted(vFaces, grid, pElem);
805 0 : vSubsetIndex.resize(vFaces.size());
806 0 : for(size_t i = 0; i < vFaces.size(); ++i)
807 0 : vSubsetIndex[i] = ish->get_subset_index(vFaces[i]);
808 0 : }
809 :
810 : try{
811 : const DimReferenceElement<dim>& rRefElem
812 0 : = ReferenceElementProvider::get<dim>(m_roid);
813 :
814 0 : DimReferenceMapping<dim, worldDim>& rMapping = ReferenceMappingProvider::get<dim, worldDim>(m_roid);
815 0 : rMapping.update(vCornerCoords);
816 :
817 : const LocalShapeFunctionSet<dim>& TrialSpace =
818 0 : LocalFiniteElementProvider::get<dim>(m_roid, LFEID(LFEID::LAGRANGE, dim, 1));
819 :
820 : // loop requested subset
821 : typename std::map<int, std::vector<BF> >::iterator it;
822 0 : for (it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); ++it)
823 : {
824 : // get subset index
825 0 : const int bndIndex = (*it).first;
826 :
827 : // get vector of BF for element
828 0 : std::vector<BF>& vBF = (*it).second;
829 :
830 : // clear vector
831 : vBF.clear();
832 :
833 : // current number of bf
834 : size_t curr_bf = 0;
835 :
836 : // loop sides of element
837 0 : for(size_t side = 0; side < vSubsetIndex.size(); ++side)
838 : {
839 : // skip non boundary sides
840 0 : if(vSubsetIndex[side] != bndIndex) continue;
841 :
842 : // number of corners of side
843 0 : const int coOfSide = rRefElem.num(dim-1, side, 0);
844 :
845 : // resize vector
846 0 : vBF.resize(curr_bf + coOfSide);
847 :
848 : // loop corners
849 0 : for(int co = 0; co < coOfSide; ++co)
850 : {
851 : // get current bf
852 : BF& bf = vBF[curr_bf];
853 :
854 : // set node id == scv this bf belongs to
855 0 : bf.nodeId = rRefElem.id(dim-1, side, 0, co);
856 :
857 : // Compute MidID for BF
858 0 : ComputeBFMidID(rRefElem, side, bf.vMidID, co);
859 :
860 : // copy corners of bf
861 0 : CopyCornerByMidID<dim, maxMid>(bf.vLocPos, bf.vMidID, m_vvLocMid, BF::numCo);
862 0 : CopyCornerByMidID<worldDim, maxMid>(bf.vGloPos, bf.vMidID, m_vvGloMid, BF::numCo);
863 :
864 : // integration point
865 0 : AveragePositions(bf.localIP, bf.vLocPos, BF::numCo);
866 0 : AveragePositions(bf.globalIP, bf.vGloPos, BF::numCo);
867 :
868 : // normal on scvf
869 0 : traits::NormalOnBF(bf.Normal, bf.vGloPos, m_vvGloMid[0]);
870 :
871 : // compute volume
872 0 : bf.Vol = VecTwoNorm(bf.Normal);
873 :
874 : // compute shapes and grads
875 0 : bf.numSH = TrialSpace.num_sh();
876 0 : TrialSpace.shapes(&(bf.vShape[0]), bf.localIP);
877 0 : TrialSpace.grads(&(bf.vLocalGrad[0]), bf.localIP);
878 :
879 : // get reference mapping
880 0 : rMapping.jacobian_transposed_inverse(bf.JtInv, bf.localIP);
881 0 : bf.detj = rMapping.sqrt_gram_det(bf.localIP);
882 :
883 : // compute global gradients
884 0 : for(size_t sh = 0 ; sh < bf.num_sh(); ++sh)
885 0 : MatVecMult(bf.vGlobalGrad[sh], bf.JtInv, bf.vLocalGrad[sh]);
886 :
887 : // increase curr_bf
888 0 : ++curr_bf;
889 : }
890 : }
891 : }
892 :
893 : }
894 0 : UG_CATCH_THROW("DimFV1Geometry: update failed.");
895 0 : }
896 :
897 : ////////////////////////////////////////////////////////////////////////////////
898 : // FV1ManifoldGeometry
899 : ////////////////////////////////////////////////////////////////////////////////
900 :
901 : template <typename TElem, int TWorldDim>
902 0 : FV1ManifoldGeometry<TElem, TWorldDim>::
903 0 : FV1ManifoldGeometry() : m_pElem(NULL), m_rRefElem(Provider<ref_elem_type>::get())
904 : {
905 : // set corners of element as local centers of nodes
906 0 : for (size_t i = 0; i < m_rRefElem.num(0); ++i)
907 : m_locMid[0][i] = m_rRefElem.corner(i);
908 :
909 : // compute local midpoints for all geometric objects with 0 < d <= dim
910 0 : for (int d = 1; d <= dim; ++d)
911 : {
912 : // loop geometric objects of dimension d
913 0 : for(size_t i = 0; i < m_rRefElem.num(d); ++i)
914 : {
915 : // set first node
916 0 : const size_t coID0 = m_rRefElem.id(d, i, 0, 0);
917 : m_locMid[d][i] = m_locMid[0][coID0];
918 :
919 : // add corner coordinates of the corners of the geometric object
920 0 : for(size_t j = 1; j < m_rRefElem.num(d, i, 0); ++j)
921 : {
922 0 : const size_t coID = m_rRefElem.id(d, i, 0, j);
923 0 : m_locMid[d][i] += m_locMid[0][coID];
924 : }
925 :
926 : // scale for correct averaging
927 0 : m_locMid[d][i] *= 1./(m_rRefElem.num(d, i, 0));
928 : }
929 : }
930 :
931 : // set up local information for Boundary Faces (bf)
932 : // each bf is associated to one corner of the element
933 0 : for (size_t i = 0; i < num_bf(); ++i)
934 : {
935 0 : m_vBF[i].nodeId = i;
936 :
937 : if (dim == 1) // RegularEdge
938 : {
939 0 : m_vBF[i].midId[0] = MidID(0, i); // set node as corner of bf
940 0 : m_vBF[i].midId[1] = MidID(dim, 0); // center of bnd element
941 :
942 : // copy local corners of bf
943 : copy_local_corners(m_vBF[i]);
944 : }
945 : else if (dim == 2) // Quadrilateral
946 : {
947 0 : m_vBF[i].midId[0] = MidID(0, i); // set node as corner of bf
948 0 : m_vBF[i].midId[1] = MidID(1, m_rRefElem.id(0, i, 1, 0)); // edge 1
949 0 : m_vBF[i].midId[2] = MidID(dim, 0); // center of bnd element
950 0 : m_vBF[i].midId[3] = MidID(1, m_rRefElem.id(0, i, 1, 1)); // edge 2
951 :
952 : // copy local corners of bf
953 : copy_local_corners(m_vBF[i]);
954 : }
955 : else {UG_ASSERT(0, "Dimension higher than 2 not implemented.");}
956 : }
957 :
958 : /////////////
959 : // Shapes
960 : /////////////
961 : // Shapes are calculated using the lower-dimensional shape functions
962 : // on the manifold element. This is the same (for Lagrange shape functions)
963 : // as using the actual shapes of the full-dimensional element.
964 0 : for (size_t i = 0; i < num_bf(); ++i)
965 : {
966 : const LocalShapeFunctionSet<ref_elem_type::dim>& TrialSpace =
967 : LocalFiniteElementProvider::get<ref_elem_type::dim>
968 : (
969 : ref_elem_type::REFERENCE_OBJECT_ID,
970 0 : LFEID(LFEID::LAGRANGE, ref_elem_type::dim, 1)
971 : );
972 :
973 : const size_t num_sh = m_numBF;
974 0 : m_vBF[i].vShape.resize(num_sh);
975 :
976 0 : TrialSpace.shapes(&(m_vBF[i].vShape[0]), m_vBF[i].local_ip());
977 : }
978 :
979 : ///////////////////////////
980 : // Copy ip pos in list
981 : ///////////////////////////
982 :
983 : // loop Boundary Faces (BF)
984 : m_vLocBFIP.clear();
985 0 : for (size_t i = 0; i < num_bf(); ++i)
986 : {
987 : // get current BF
988 : const BF& rBF = bf(i);
989 :
990 : // loop ips
991 0 : m_vLocBFIP.push_back(rBF.local_ip());
992 : }
993 0 : }
994 :
995 :
996 : /// update data for given element
997 : template <typename TElem, int TWorldDim>
998 0 : void FV1ManifoldGeometry<TElem, TWorldDim>::
999 : update(GridObject* elem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
1000 : {
1001 : // if already update for this element, do nothing
1002 0 : if (m_pElem == elem) return;
1003 0 : else m_pElem = elem;
1004 :
1005 : // remember global position of nodes
1006 0 : for (size_t i = 0; i < m_rRefElem.num(0); ++i)
1007 0 : m_gloMid[0][i] = vCornerCoords[i];
1008 :
1009 : // compute global midpoints for all the other geometric objects (with 0 < d <= dim)
1010 0 : for (int d = 1; d <= dim; ++d)
1011 : {
1012 : // loop geometric objects of dimension d
1013 0 : for (size_t i = 0; i < m_rRefElem.num(d); ++i)
1014 : {
1015 : // set first node
1016 0 : const size_t coID0 = m_rRefElem.id(d, i, 0, 0);
1017 : m_gloMid[d][i] = m_gloMid[0][coID0];
1018 :
1019 : // add corner coordinates of the corners of the geometric object
1020 0 : for (size_t j = 1; j < m_rRefElem.num(d, i, 0); ++j)
1021 : {
1022 0 : const size_t coID = m_rRefElem.id(d, i, 0, j);
1023 : m_gloMid[d][i] += m_gloMid[0][coID];
1024 : }
1025 :
1026 : // scale for correct averaging
1027 0 : m_gloMid[d][i] *= 1./(m_rRefElem.num(d, i, 0));
1028 : }
1029 : }
1030 :
1031 : // set global integration points
1032 0 : for (size_t i = 0; i < num_bf(); ++i)
1033 : {
1034 : // copy global corners of bf
1035 0 : copy_global_corners(m_vBF[i]);
1036 : }
1037 :
1038 : // compute size of BFs
1039 0 : for (size_t i = 0; i < num_bf(); ++i)
1040 : {
1041 : // copy global corners
1042 0 : copy_global_corners(m_vBF[i]);
1043 :
1044 : // compute volume of bf
1045 0 : m_vBF[i].vol = ElementSize<bf_type, worldDim>(m_vBF[i].vGloPos);
1046 : }
1047 :
1048 : ///////////////////////////
1049 : // Copy ip pos in list
1050 : ///////////////////////////
1051 :
1052 : // loop Boundary Faces (BF)
1053 : m_vGlobBFIP.clear();
1054 0 : for (size_t i = 0; i < num_bf(); ++i)
1055 : {
1056 : // get current BF
1057 : const BF& rBF = bf(i);
1058 :
1059 0 : m_vGlobBFIP.push_back(rBF.global_ip());
1060 : }
1061 : }
1062 :
1063 : //////////////////////
1064 : // FV1Geometry
1065 :
1066 : template class FV1Geometry_gen<RegularEdge, 1, false>;
1067 : template class FV1Geometry_gen<RegularEdge, 2, false>;
1068 : template class FV1Geometry_gen<RegularEdge, 3, false>;
1069 :
1070 : template class FV1Geometry_gen<Triangle, 2, false>;
1071 : template class FV1Geometry_gen<Triangle, 3, false>;
1072 :
1073 : template class FV1Geometry_gen<Quadrilateral, 2, false>;
1074 : template class FV1Geometry_gen<Quadrilateral, 3, false>;
1075 :
1076 : template class FV1Geometry_gen<Tetrahedron, 3, false>;
1077 : template class FV1Geometry_gen<Prism, 3, false>;
1078 : template class FV1Geometry_gen<Pyramid, 3, false>;
1079 : template class FV1Geometry_gen<Hexahedron, 3, false>;
1080 : template class FV1Geometry_gen<Octahedron, 3, false>;
1081 :
1082 : //////////////////////
1083 : // 'condensed' FV1Geometry
1084 :
1085 : template class FV1Geometry_gen<RegularEdge, 1, true>;
1086 : template class FV1Geometry_gen<RegularEdge, 2, true>;
1087 : template class FV1Geometry_gen<RegularEdge, 3, true>;
1088 :
1089 : template class FV1Geometry_gen<Triangle, 2, true>;
1090 : template class FV1Geometry_gen<Triangle, 3, true>;
1091 :
1092 : template class FV1Geometry_gen<Quadrilateral, 2, true>;
1093 : template class FV1Geometry_gen<Quadrilateral, 3, true>;
1094 :
1095 : template class FV1Geometry_gen<Tetrahedron, 3, true>;
1096 : template class FV1Geometry_gen<Prism, 3, true>;
1097 : template class FV1Geometry_gen<Pyramid, 3, true>;
1098 : template class FV1Geometry_gen<Hexahedron, 3, true>;
1099 : template class FV1Geometry_gen<Octahedron, 3, true>;
1100 :
1101 : //////////////////////
1102 : // DimFV1Geometry
1103 : template class DimFV1Geometry<1, 1>;
1104 : template class DimFV1Geometry<1, 2>;
1105 : template class DimFV1Geometry<1, 3>;
1106 :
1107 : template class DimFV1Geometry<2, 2>;
1108 : template class DimFV1Geometry<2, 3>;
1109 :
1110 : template class DimFV1Geometry<3, 3>;
1111 :
1112 : //////////////////////
1113 : // Manifold
1114 : template class FV1ManifoldGeometry<RegularEdge, 2>;
1115 : template class FV1ManifoldGeometry<Triangle, 3>;
1116 : template class FV1ManifoldGeometry<Quadrilateral, 3>;
1117 :
1118 : } // end namespace ug
|
