Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Christian Wehner
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : /*
34 : * Node centered finite volume geometry for Crouzeix-Raviart-Elements
35 : */
36 :
37 : #include "common/util/provider.h"
38 : #include "fvcr_geom.h"
39 : #include "lib_disc/reference_element/reference_element.h"
40 : #include "lib_disc/reference_element/reference_mapping.h"
41 : #include "lib_disc/reference_element/reference_mapping_provider.h"
42 : #include "lib_disc/quadrature/quadrature.h"
43 :
44 : namespace ug{
45 :
46 : ////////////////////////////////////////////////////////////////////////////////
47 : ////////////////////////////////////////////////////////////////////////////////
48 : // Dim-dependent Crouzeix-Raviart Finite Volume Geometry
49 : ////////////////////////////////////////////////////////////////////////////////
50 : ////////////////////////////////////////////////////////////////////////////////
51 :
52 : template <int TDim, int TWorldDim>
53 0 : void DimCRFVGeometry<TDim, TWorldDim>::
54 : update_local_data()
55 : {
56 : // get reference element
57 : try{
58 : m_rRefElem
59 0 : = ReferenceElementProvider::get<dim>(m_roid);
60 :
61 : // set number of scvf / scv of this roid
62 0 : m_numSCV = m_rRefElem.num(dim-1); // number of faces
63 0 : m_numSCVF = m_rRefElem.num(1); // number of edges
64 :
65 : // compute barycenter coordinates
66 : localBary = m_rRefElem.corner(0);
67 0 : for (size_t j=1;j<m_rRefElem.num(0);j++){
68 : localBary+=m_rRefElem.corner(j);
69 : }
70 0 : localBary*=1./(number)m_rRefElem.num(0);
71 :
72 : // set up local informations for SubControlVolumeFaces (scvf)
73 : // each scvf is associated to one vertex (2d) / edge (3d) of the element
74 0 : for(size_t i = 0; i < m_numSCVF; ++i)
75 : {
76 : // this scvf separates the given edges/faces
77 0 : m_vSCVF[i].From = m_rRefElem.id(dim-2, i, dim-1, 0);// todo handle dim==1
78 0 : m_vSCVF[i].To = m_rRefElem.id(dim-2, i, dim-1, 1);
79 :
80 0 : for (size_t j=0;j<m_vSCVF[i].numCo-1;j++){
81 0 : m_vSCVF[i].vLocPos[j]=m_rRefElem.corner(m_rRefElem.id(dim-2,i,0,j));
82 : }
83 :
84 : m_vSCVF[i].vLocPos[m_vSCVF[i].numCo-1]=localBary;
85 0 : AveragePositions(m_vSCVF[i].localIP, m_vSCVF[i].vLocPos, m_vSCVF[i].numCo);
86 : }
87 :
88 : // set up local informations for SubControlVolumes (scv)
89 : // each scv is associated to one edge(2d) / face(3d) of the element
90 0 : for(size_t i = 0; i < m_numSCV; ++i)
91 : {
92 : // store associated node
93 0 : m_vSCV[i].nodeID = i;
94 :
95 0 : m_vSCV[i].numCorners = m_rRefElem.num(dim-1,i,0)+1;
96 0 : for (int j=0;j<m_vSCV[i].numCorners-1;j++){
97 0 : m_vSCV[i].vLocPos[m_vSCV[i].numCorners-2-j]=m_rRefElem.corner(m_rRefElem.id(dim-1,i,0,j));
98 : }
99 0 : AveragePositions(m_vLocUnkCoords[i], m_vSCV[i].vLocPos, m_vSCV[i].numCorners-1);
100 : m_vSCV[i].vLocIP = m_vLocUnkCoords[i];
101 : m_vSCV[i].vLocPos[m_vSCV[i].numCorners-1]=localBary;
102 : }
103 :
104 : /////////////////////////
105 : // Shapes and Derivatives
106 : /////////////////////////
107 :
108 : const LocalShapeFunctionSet<dim>& rTrialSpace =
109 0 : LocalFiniteElementProvider::get<dim>(m_roid, LFEID(LFEID::CROUZEIX_RAVIART, dim, 1));
110 :
111 0 : m_nsh = rTrialSpace.num_sh();
112 :
113 0 : for(size_t i = 0; i < m_numSCVF; ++i)
114 : {
115 0 : m_vSCVF[i].numSH = rTrialSpace.num_sh();
116 0 : rTrialSpace.shapes(&(m_vSCVF[i].vShape[0]), m_vSCVF[i].local_ip());
117 0 : rTrialSpace.grads(&(m_vSCVF[i].vLocalGrad[0]), m_vSCVF[i].local_ip());
118 : }
119 :
120 0 : for(size_t i = 0; i < m_numSCV; ++i)
121 : {
122 0 : m_vSCV[i].numSH = rTrialSpace.num_sh();
123 0 : rTrialSpace.shapes(&(m_vSCV[i].vShape[0]), m_vSCV[i].local_ip());
124 0 : rTrialSpace.grads(&(m_vSCV[i].vLocalGrad[0]), m_vSCV[i].local_ip());
125 : }
126 :
127 : }
128 0 : UG_CATCH_THROW("DimCRFVGeometry: update failed.");
129 :
130 : // copy ip positions in a list for Sub Control Volumes Faces (SCVF)
131 0 : for(size_t i = 0; i < m_numSCVF; ++i)
132 : m_vLocSCVF_IP[i] = scvf(i).local_ip();
133 :
134 0 : m_numConstrainedDofs = 0;
135 0 : m_numConstrainedSCVF = 0;
136 0 : }
137 :
138 :
139 : /// update data for given element
140 : template <int TDim, int TWorldDim>
141 0 : void DimCRFVGeometry<TDim, TWorldDim>::
142 : update(GridObject* pElem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
143 : {
144 : // If already update for this element, do nothing
145 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
146 :
147 : // refresh local data, if different roid given
148 0 : if(m_roid != pElem->reference_object_id())
149 : {
150 : // remember new roid
151 0 : m_roid = (ReferenceObjectID) pElem->reference_object_id();
152 :
153 : // update local data
154 0 : update_local_data();
155 : }
156 :
157 : // get reference element
158 : try{
159 : const DimReferenceElement<dim>& m_rRefElem
160 0 : = ReferenceElementProvider::get<dim>(m_roid);
161 :
162 : // compute barycenter coordinates
163 : globalBary = vCornerCoords[0];
164 0 : for (size_t j=1;j<m_rRefElem.num(0);j++){
165 0 : globalBary+=vCornerCoords[j];
166 : }
167 0 : globalBary*=1./(number)m_rRefElem.num(0);
168 :
169 : // compute global informations for scvf
170 0 : for(size_t i = 0; i < num_scvf(); ++i)
171 : {
172 0 : for (size_t j=0;j<m_vSCVF[i].numCo-1;j++){
173 0 : m_vSCVF[i].vGloPos[j]=vCornerCoords[m_rRefElem.id(dim-2,i,0,j)];
174 : }
175 : m_vSCVF[i].vGloPos[m_vSCVF[i].numCo-1]=globalBary;
176 0 : AveragePositions(m_vSCVF[i].globalIP, m_vSCVF[i].vGloPos, m_vSCVF[i].numCo);
177 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[i].Normal,m_vSCVF[i].vGloPos);// face_type0 identical to scvf type
178 : }
179 :
180 : // compute size of scv
181 0 : for(size_t i = 0; i < num_scv(); ++i)
182 : {
183 : // side nodes in reverse order to fulfill standard element order
184 0 : for (int j=0;j<m_vSCV[i].numCorners-1;j++){
185 0 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-2-j]=vCornerCoords[m_rRefElem.id(dim-1,i,0,j)];
186 : }
187 0 : AveragePositions(m_vGlobUnkCoords[i], m_vSCV[i].vGloPos, m_vSCV[i].numCorners-1);
188 : m_vSCV[i].vGlobIP = m_vGlobUnkCoords[i];
189 :
190 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-1]=globalBary;
191 : // compute volume of scv and normal to associated element face
192 : //CRSCVSizeAndNormal<dim>(m_vSCV[i].Vol,m_vSCV[i].Normal,m_vSCV[i].vGloPos,m_vSCV[i].numCorners);
193 0 : if (m_vSCV[i].numCorners-1==dim){
194 0 : m_vSCV[i].Vol = ElementSize<scv_type0,worldDim>(m_vSCV[i].vGloPos);
195 0 : ElementNormal<face_type0, worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
196 : } else { // m_vSCV[i].numCorners-2==dim , only possible in 3d (pyramid)
197 0 : m_vSCV[i].Vol = ElementSize<scv_type1,worldDim>(m_vSCV[i].vGloPos);
198 0 : ElementNormal<face_type1,worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
199 : };
200 : // nodes are in reverse order therefore reverse sign to get outward normal
201 : m_vSCV[i].Normal*=-1;
202 : }
203 :
204 : // get reference mapping
205 0 : DimReferenceMapping<dim, worldDim>& rMapping = ReferenceMappingProvider::get<dim, worldDim>(m_roid);
206 0 : rMapping.update(vCornerCoords);
207 :
208 : //\todo compute with on virt. call
209 : // compute jacobian for linear mapping
210 0 : if(rMapping.is_linear())
211 : {
212 : MathMatrix<worldDim,dim> JtInv;
213 0 : rMapping.jacobian_transposed_inverse(JtInv, m_vSCVF[0].local_ip());
214 0 : const number detJ = rMapping.sqrt_gram_det(m_vSCVF[0].local_ip());
215 :
216 0 : for(size_t i = 0; i < num_scvf(); ++i)
217 : {
218 : m_vSCVF[i].JtInv = JtInv;
219 0 : m_vSCVF[i].detj = detJ;
220 : }
221 :
222 0 : for(size_t i = 0; i < num_scv(); ++i)
223 : {
224 : m_vSCV[i].JtInv = JtInv;
225 0 : m_vSCV[i].detj = detJ;
226 : }
227 : }
228 : // else compute jacobian for each integration point
229 : else
230 : {
231 0 : for(size_t i = 0; i < num_scvf(); ++i)
232 : {
233 0 : rMapping.jacobian_transposed_inverse(m_vSCVF[i].JtInv, m_vSCVF[i].local_ip());
234 0 : m_vSCVF[i].detj = rMapping.sqrt_gram_det(m_vSCVF[i].local_ip());
235 : }
236 0 : for(size_t i = 0; i < num_scv(); ++i)
237 : {
238 0 : rMapping.jacobian_transposed_inverse(m_vSCV[i].JtInv, m_vSCV[i].local_ip());
239 0 : m_vSCV[i].detj = rMapping.sqrt_gram_det(m_vSCV[i].local_ip());
240 : }
241 : }
242 :
243 : // compute global gradients
244 0 : for(size_t i = 0; i < num_scvf(); ++i)
245 0 : for(size_t sh = 0; sh < scvf(i).num_sh(); ++sh)
246 : MatVecMult(m_vSCVF[i].vGlobalGrad[sh], m_vSCVF[i].JtInv, m_vSCVF[i].vLocalGrad[sh]);
247 :
248 0 : for(size_t i = 0; i < num_scv(); ++i)
249 0 : for(size_t sh = 0; sh < scv(i).num_sh(); ++sh)
250 : MatVecMult(m_vSCV[i].vGlobalGrad[sh], m_vSCV[i].JtInv, m_vSCV[i].vLocalGrad[sh]);
251 :
252 : // copy ip points in list (SCVF)
253 0 : for(size_t i = 0; i < num_scvf(); ++i)
254 : m_vGlobSCVF_IP[i] = scvf(i).global_ip();
255 :
256 : }
257 0 : UG_CATCH_THROW("DimCRFVGeometry: update failed.");
258 :
259 : // if no boundary subsets required, return
260 0 : if(num_boundary_subsets() == 0 || ish == NULL) return;
261 0 : else update_boundary_faces(pElem, vCornerCoords, ish);
262 : }
263 :
264 : /// update data checking for hanging nodes for given element
265 : template <int TDim, int TWorldDim>
266 0 : void DimCRFVGeometry<TDim, TWorldDim>::
267 : update_hanging(GridObject* pElem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish,bool keepSCV,bool keepSCVF)
268 : {
269 : // If already update for this element, do nothing
270 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
271 :
272 : // get grid
273 0 : Grid& grid = *(ish->grid());
274 :
275 : // refresh local data, if different roid given
276 0 : if(m_roid != pElem->reference_object_id())
277 : {
278 : // remember new roid
279 0 : m_roid = (ReferenceObjectID) pElem->reference_object_id();
280 :
281 : // update local data
282 0 : update_local_data();
283 :
284 0 : m_numDofs = num_scv();
285 : }
286 :
287 : const LocalShapeFunctionSet<dim>& rTrialSpace =
288 0 : LocalFiniteElementProvider::get<dim>(m_roid, LFEID(LFEID::CROUZEIX_RAVIART, dim, 1));
289 :
290 : // compute barycenter coordinates
291 : globalBary = vCornerCoords[0];
292 0 : for (size_t j=1;j<m_rRefElem.num(0);j++){
293 0 : globalBary+=vCornerCoords[j];
294 : }
295 0 : globalBary*=1./(number)m_rRefElem.num(0);
296 :
297 : // compute global informations for scvf
298 0 : for(size_t i = 0; i < num_scvf(); ++i)
299 : {
300 0 : for (size_t j=0;j<m_vSCVF[i].numCo-1;j++){
301 0 : m_vSCVF[i].vGloPos[j]=vCornerCoords[m_rRefElem.id(dim-2,i,0,j)];
302 : }
303 : m_vSCVF[i].vGloPos[m_vSCVF[i].numCo-1]=globalBary;
304 0 : AveragePositions(m_vSCVF[i].globalIP, m_vSCVF[i].vGloPos, m_vSCVF[i].numCo);
305 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[i].Normal,m_vSCVF[i].vGloPos);// face_type0 identical to scvf type
306 : }
307 :
308 : // compute size of scv
309 0 : for(size_t i = 0; i < num_scv(); ++i)
310 : {
311 : // side nodes in reverse order to fulfill standard element order
312 0 : for (int j=0;j<m_vSCV[i].numCorners-1;j++){
313 0 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-2-j]=vCornerCoords[m_rRefElem.id(dim-1,i,0,j)];
314 : }
315 0 : AveragePositions(m_vGlobUnkCoords[i], m_vSCV[i].vGloPos, m_vSCV[i].numCorners-1);
316 : m_vSCV[i].vGlobIP = m_vGlobUnkCoords[i];
317 :
318 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-1]=globalBary;
319 : // compute volume of scv and normal to associated element face
320 : //CRSCVSizeAndNormal<dim>(m_vSCV[i].Vol,m_vSCV[i].Normal,m_vSCV[i].vGloPos,m_vSCV[i].numCorners);
321 0 : if (m_vSCV[i].numCorners-1==dim){
322 0 : m_vSCV[i].Vol = ElementSize<scv_type0,worldDim>(m_vSCV[i].vGloPos);
323 0 : ElementNormal<face_type0, worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
324 : } else { // m_vSCV[i].numCorners-2==dim , only possible in 3d (pyramid)
325 0 : m_vSCV[i].Vol = ElementSize<scv_type1,worldDim>(m_vSCV[i].vGloPos);
326 0 : ElementNormal<face_type1,worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
327 : };
328 : // nodes are in reverse order therefore reverse sign to get outward normal
329 : m_vSCV[i].Normal*=-1;
330 : }
331 :
332 : // get reference mapping
333 0 : DimReferenceMapping<dim, worldDim>& rMapping = ReferenceMappingProvider::get<dim, worldDim>(m_roid);
334 0 : rMapping.update(vCornerCoords);
335 :
336 : // check for hanging nodes
337 : if (dim==2){
338 : std::vector<Edge*> vEdges;
339 0 : CollectEdgesSorted(vEdges, grid, pElem);
340 0 : for(size_t side = 0; side < vEdges.size(); ++side){
341 0 : ConstrainingEdge* constrainingObj = dynamic_cast<ConstrainingEdge*>(vEdges[side]);
342 0 : if(constrainingObj == NULL) continue;
343 :
344 : // found constraining edge
345 : MathVector<worldDim> globalMidpoint = m_vSCV[side].vGlobIP;
346 : MathVector<dim> localMidpoint = m_vSCV[side].vLocIP;
347 : // get edge corners
348 : size_t edgeCo[2];
349 0 : for (size_t j=0;j<2;j++) edgeCo[j] = m_rRefElem.id(1,side,0,j);
350 : // compute dof positions on constraining edge,
351 : // replace dof "side" with first and insert second at the end
352 : // set up new scvs
353 : // keepSCV parameter specifies if scv "side" gets replaced by new one
354 : size_t ind=side;
355 : size_t keepOffset=0;
356 : MathVector<worldDim> normal = m_vSCV[side].Normal;
357 : normal*=0.5;
358 0 : number vol = 0.5*m_vSCV[side].Vol;
359 0 : if (keepSCV){
360 0 : ind=m_numSCV+1;
361 : keepOffset=1;
362 : }
363 0 : for (int d=0;d<worldDim;d++)
364 0 : m_vGlobUnkCoords[ind][d] = 0.5 * (vCornerCoords[edgeCo[0]][d] + globalMidpoint[d]);
365 0 : for (int d=0;d<dim;d++)
366 0 : m_vLocUnkCoords[ind][d] = 0.5 * (m_rRefElem.corner(edgeCo[0])[d] + localMidpoint[d]);
367 0 : for (int d=0;d<worldDim;d++)
368 0 : m_vGlobUnkCoords[m_numSCV][d] = 0.5 * (vCornerCoords[edgeCo[1]][d] + globalMidpoint[d]);
369 0 : for (int d=0;d<dim;d++)
370 0 : m_vLocUnkCoords[m_numSCV][d] = 0.5 * (m_rRefElem.corner(edgeCo[1])[d] + localMidpoint[d]);
371 : // handle corresponding scvfs
372 : // if keepSCVF copy them into constrained scvf array
373 0 : if (keepSCVF){
374 0 : for (size_t j=0;j<2;j++){
375 0 : m_vConstrainedSCVF[m_numConstrainedSCVF+j] = m_vSCVF[edgeCo[j]];
376 0 : if (m_vConstrainedSCVF[edgeCo[j]].To==side){
377 0 : m_vConstrainedSCVF[edgeCo[j]].From=side;
378 : m_vConstrainedSCVF[edgeCo[j]].Normal*=-1;
379 : }
380 : }
381 0 : m_numConstrainedSCVF += 2;
382 : }
383 0 : for (size_t j=0;j<2;j++){
384 0 : if (m_vSCVF[edgeCo[j]].From==side) m_vSCVF[edgeCo[j]].From=m_numDofs+j;
385 0 : if (m_vSCVF[edgeCo[j]].To==side) m_vSCVF[edgeCo[j]].To=m_numDofs+j;
386 : }
387 : m_vSCV[ind].Normal = normal;
388 0 : m_vSCV[ind].Vol = vol;
389 0 : m_vSCV[ind].nodeID = m_numDofs;
390 : m_vSCV[ind].vGlobIP = m_vGlobUnkCoords[ind];
391 : m_vSCV[ind].vLocIP = m_vLocUnkCoords[ind];
392 0 : m_vSCV[ind].numSH = rTrialSpace.num_sh();
393 0 : m_vSCV[ind].vGloPos[0]=vCornerCoords[edgeCo[0]];
394 : m_vSCV[ind].vGloPos[1]=globalMidpoint;
395 : m_vSCV[ind].vGloPos[2]=globalBary;
396 : m_vSCV[ind].vLocPos[0]=m_rRefElem.corner(edgeCo[0]);
397 : m_vSCV[ind].vLocPos[1]=localMidpoint;
398 : m_vSCV[ind].vLocPos[2]=localBary;
399 0 : m_vSCV[ind].numCorners = 3;
400 0 : rTrialSpace.shapes(&(m_vSCV[ind].vShape[0]), m_vSCV[ind].local_ip());
401 0 : rTrialSpace.grads(&(m_vSCV[ind].vLocalGrad[0]), m_vSCV[ind].local_ip());
402 : // second scv inserted at the end
403 0 : m_vSCV[m_numSCV].Normal = normal;
404 0 : m_vSCV[m_numSCV].Vol = vol;
405 0 : m_vSCV[m_numSCV].nodeID = m_numDofs+1;
406 : m_vSCV[m_numSCV].vGlobIP = m_vGlobUnkCoords[m_numSCV];
407 : m_vSCV[m_numSCV].vLocIP = m_vLocUnkCoords[m_numSCV];
408 0 : m_vSCV[m_numSCV].numSH = rTrialSpace.num_sh();
409 0 : m_vSCV[m_numSCV].vGloPos[0]=vCornerCoords[edgeCo[1]];
410 : m_vSCV[m_numSCV].vGloPos[1]=globalMidpoint;
411 : m_vSCV[m_numSCV].vGloPos[2]=globalBary;
412 : m_vSCV[m_numSCV].vLocPos[0]=m_rRefElem.corner(edgeCo[1]);
413 : m_vSCV[m_numSCV].vLocPos[1]=localMidpoint;
414 : m_vSCV[m_numSCV].vLocPos[2]=localBary;
415 0 : m_vSCV[m_numSCV].numCorners = 3;
416 0 : rTrialSpace.shapes(&(m_vSCV[m_numSCV].vShape[0]), m_vSCV[m_numSCV].local_ip());
417 0 : rTrialSpace.grads(&(m_vSCV[m_numSCV].vLocalGrad[0]), m_vSCV[m_numSCV].local_ip());
418 : // insert new scvf
419 0 : m_vSCVF[m_numSCVF].From = m_numDofs;
420 0 : m_vSCVF[m_numSCVF].To = m_numDofs+1;
421 : m_vSCVF[m_numSCVF].vLocPos[0] = localMidpoint;
422 : m_vSCVF[m_numSCVF].vLocPos[1] = localBary;
423 : m_vSCVF[m_numSCVF].vGloPos[0] = globalMidpoint;
424 : m_vSCVF[m_numSCVF].vGloPos[1] = globalBary;
425 0 : m_vSCVF[m_numSCVF].numSH = rTrialSpace.num_sh();
426 0 : for (int d=0;d<dim;d++) m_vSCVF[m_numSCVF].localIP[d] = 0.5*(localMidpoint[d] + localBary[d]);
427 0 : for (int d=0;d<worldDim;d++) m_vSCVF[m_numSCVF].globalIP[d] = 0.5*(globalMidpoint[d] + globalBary[d]);
428 0 : rTrialSpace.shapes(&(m_vSCVF[m_numSCVF].vShape[0]), m_vSCVF[m_numSCVF].local_ip());
429 0 : rTrialSpace.grads(&(m_vSCVF[m_numSCVF].vLocalGrad[0]), m_vSCVF[m_numSCVF].local_ip());
430 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[m_numSCVF].Normal,m_vSCVF[m_numSCVF].vGloPos);
431 : // insert new constrained dof object
432 0 : m_vCD[m_numConstrainedDofs].i = side;
433 0 : m_vCD[m_numConstrainedDofs].numConstrainingDofs = 2;
434 0 : m_vCD[m_numConstrainedDofs].cDofInd[0] = m_numDofs;
435 0 : m_vCD[m_numConstrainedDofs].cDofInd[1] = m_numDofs+1;
436 0 : m_vCD[m_numConstrainedDofs].cDofWeights[0] = 0.5;
437 0 : m_vCD[m_numConstrainedDofs].cDofWeights[1] = 0.5;
438 0 : m_numSCV+=1+keepOffset;
439 0 : m_numSCVF+=1;
440 0 : m_numDofs+=2;
441 0 : m_numConstrainedDofs+=1;
442 0 : m_roid = ROID_UNKNOWN;
443 : }
444 0 : } else {
445 : // dim == 3
446 : std::vector<Face*> vFaces;
447 0 : CollectFacesSorted(vFaces, grid, pElem);
448 : handledEdges.clear();
449 0 : for(size_t face = 0; face < vFaces.size(); ++face){
450 0 : ConstrainingFace* constrainingObj = dynamic_cast<ConstrainingFace*>(vFaces[face]);
451 0 : if(constrainingObj == NULL) continue;
452 : // found constraining face
453 : MathVector<worldDim> globalMidpoint = m_vSCV[face].vGlobIP;
454 : MathVector<dim> localMidpoint = m_vSCV[face].vLocIP;
455 0 : number faceVol = m_vSCV[face].Vol;
456 : MathVector<worldDim> faceNormal = m_vSCV[face].Normal;
457 : // get face corners and edges
458 : size_t faceCo[4];
459 : size_t faceEdge[4];
460 : // number of corners of face = number of edges
461 : size_t numFaceCo = m_rRefElem.num(2,face,0);
462 0 : for (size_t j=0;j<numFaceCo;j++) faceCo[j] = m_rRefElem.id(2,face,0,j);
463 0 : for (size_t j=0;j<numFaceCo;j++) faceEdge[j] = m_rRefElem.id(2,face,1,j);
464 : number volSum=0;
465 : size_t keepOffset=0;
466 0 : if (keepSCV) keepOffset=1;
467 : // compute coordinates of each face and fill scv values
468 0 : for (size_t i=0;i<numFaceCo;i++){
469 0 : size_t co = faceCo[i];
470 : size_t nOfEdges=0;
471 : size_t nbEdges[2];
472 : // find 2 edges in face belonging to node
473 0 : for (size_t j=0;j<m_rRefElem.num(0,co,1);j++){
474 0 : size_t candidate = m_rRefElem.id(0,co,1,j);
475 : bool found = false;
476 0 : for (size_t k=0;k<numFaceCo;k++){
477 0 : if (faceEdge[k]==candidate){
478 : found = true;
479 : break;
480 : }
481 : }
482 0 : if (found==true){
483 0 : nbEdges[nOfEdges] = candidate;
484 0 : nOfEdges++;
485 0 : if (nOfEdges==2) break;
486 : }
487 : }
488 : // in triangular case switch edges if necessary for correct orientation
489 0 : if (numFaceCo==3){
490 0 : if (faceEdge[i]==nbEdges[1]){
491 0 : nbEdges[1] = nbEdges[0];
492 0 : nbEdges[0] = faceEdge[i];
493 : }
494 : }
495 : // keepSCV parameter specifies if scv "side" gets replaced by new one
496 0 : size_t ind = m_numSCV+i-1+keepOffset;
497 0 : if (i==0){
498 0 : if (keepSCV) ind=m_numSCV;
499 : else ind = face;
500 : }
501 : // others are inserted at the end
502 0 : m_vSCV[ind].vGloPos[0] = vCornerCoords[co];
503 : m_vSCV[ind].vLocPos[0] = m_rRefElem.corner(co);
504 0 : for (int d=0;d<worldDim;d++){
505 : // edge 0 midpoint
506 0 : m_vSCV[ind].vGloPos[1][d] = 0.5 * ( vCornerCoords[m_rRefElem.id(1,nbEdges[0],0,0)][d] + vCornerCoords[m_rRefElem.id(1,nbEdges[0],0,1)][d] );
507 0 : m_vSCV[ind].vLocPos[1][d] = 0.5 * ( m_rRefElem.corner(m_rRefElem.id(1,nbEdges[0],0,0))[d] + m_rRefElem.corner(m_rRefElem.id(1,nbEdges[0],0,1))[d] );
508 : // edge 1 midpoint
509 0 : m_vSCV[ind].vGloPos[numFaceCo-1][d] = 0.5 * ( vCornerCoords[m_rRefElem.id(1,nbEdges[1],0,0)][d] + vCornerCoords[m_rRefElem.id(1,nbEdges[1],0,1)][d] );
510 0 : m_vSCV[ind].vLocPos[numFaceCo-1][d] = 0.5 * ( m_rRefElem.corner(m_rRefElem.id(1,nbEdges[1],0,0))[d] + m_rRefElem.corner(m_rRefElem.id(1,nbEdges[1],0,1))[d] );
511 : }
512 0 : if (numFaceCo==4) m_vSCV[ind].vGloPos[2] = globalMidpoint;
513 : m_vSCV[ind].vGloPos[numFaceCo] = globalBary;
514 : m_vSCV[ind].vLocPos[numFaceCo] = localBary;
515 0 : m_vSCV[ind].numCorners = numFaceCo + 1;
516 0 : AveragePositions(m_vGlobUnkCoords[ind], m_vSCV[ind].vGloPos, m_vSCV[ind].numCorners-1);
517 0 : AveragePositions(m_vLocUnkCoords[ind], m_vSCV[ind].vLocPos, m_vSCV[ind].numCorners-1);
518 : m_vSCV[ind].vLocIP = m_vLocUnkCoords[ind];
519 : m_vSCV[ind].vGlobIP = m_vGlobUnkCoords[ind];
520 0 : m_vSCV[ind].numSH = rTrialSpace.num_sh();
521 0 : if (numFaceCo==3) m_vSCV[ind].Vol = ElementSize<scv_type0,worldDim>(m_vSCV[ind].vGloPos);
522 0 : else m_vSCV[ind].Vol = ElementSize<scv_type1,worldDim>(m_vSCV[ind].vGloPos);
523 0 : if (m_vSCV[ind].Vol<0) m_vSCV[ind].Vol *= -1;
524 0 : volSum+=m_vSCV[ind].Vol;
525 : m_vSCV[ind].Normal = faceNormal;
526 0 : m_vSCV[ind].Normal *= (number) m_vSCV[ind].Vol / faceVol;
527 0 : rTrialSpace.shapes(&(m_vSCV[ind].vShape[0]), m_vSCV[ind].local_ip());
528 0 : rTrialSpace.grads(&(m_vSCV[ind].vLocalGrad[0]), m_vSCV[ind].local_ip());
529 0 : m_vSCV[ind].nodeID = m_numDofs+i;
530 : }
531 : // compute inner scv in triangular case
532 0 : if (numFaceCo==3){
533 0 : size_t ind = m_numSCV+2+keepOffset;
534 0 : m_vSCV[ind].Vol = faceVol - volSum;
535 0 : m_vSCV[ind].nodeID=m_numDofs+3;
536 : m_vSCV[ind].Normal = faceNormal;
537 0 : m_vSCV[ind].Normal *= (number) m_vSCV[ind].Vol / faceVol;
538 : m_vSCV[ind].vGlobIP = m_vSCV[face].vGloPos[1];
539 : m_vSCV[ind].vLocIP = m_vSCV[face].vLocPos[1];
540 0 : for (size_t j=0;j<2;j++){
541 0 : m_vSCV[ind].vGlobIP += m_vSCV[m_numSCV+j].vGloPos[1];
542 : m_vSCV[ind].vLocIP += m_vSCV[m_numSCV+j].vLocPos[1];
543 : }
544 : m_vSCV[ind].vGlobIP *= (number)1.0/3.0;
545 : m_vSCV[ind].vLocIP *= (number)1.0/3.0;
546 : m_vGlobUnkCoords[ind] = m_vSCV[ind].vGlobIP;
547 : m_vLocUnkCoords[ind] = m_vSCV[ind].vLocIP;
548 0 : m_vSCV[ind].numCorners = numFaceCo + 1;
549 0 : m_vSCV[ind].numSH = rTrialSpace.num_sh();
550 0 : rTrialSpace.shapes(&(m_vSCV[ind].vShape[0]), m_vSCV[ind].local_ip());
551 0 : rTrialSpace.grads(&(m_vSCV[ind].vLocalGrad[0]), m_vSCV[ind].local_ip());
552 : }
553 : // copy scvfs into constrained scvfs
554 0 : if (keepSCVF){
555 0 : for (size_t i=0;i<numFaceCo;i++){
556 0 : size_t edge = faceEdge[i];
557 0 : m_vConstrainedSCVF[m_numConstrainedSCVF+i] = m_vSCVF[edge];
558 0 : if (m_vConstrainedSCVF[m_numConstrainedSCVF+i].To==face){
559 0 : m_vConstrainedSCVF[m_numConstrainedSCVF+i].From=face;
560 : m_vConstrainedSCVF[m_numConstrainedSCVF+i].Normal*=-1;
561 : }
562 : }
563 0 : m_numConstrainedSCVF += numFaceCo;
564 : }
565 : // insert new scvfs, first the ones associated to edges of face
566 0 : for (size_t i=0;i<numFaceCo;i++){
567 0 : size_t edge = faceEdge[i];
568 0 : size_t from = m_vSCVF[edge].From;
569 0 : size_t to = m_vSCVF[edge].To;
570 : MathVector<worldDim> normal = m_vSCVF[edge].Normal;
571 : normal*=0.5;
572 : size_t edgeCo[2];
573 : size_t scvID[2];
574 0 : for (size_t j=0;j<2;j++){
575 0 : edgeCo[j] = m_rRefElem.id(1,edge,0,j);
576 : // find corresponding face vertex (= corresponding scv index)
577 0 : for (size_t k=0;k<numFaceCo;k++){
578 0 : if (faceCo[k]==edgeCo[j]){
579 0 : scvID[j] = m_numDofs+k;
580 0 : break;
581 : }
582 : }
583 : }
584 : // look if edge has already been handled
585 0 : if (m_numConstrainedDofs>0){
586 : bool found=false;
587 0 : for (size_t j=0;j<handledEdges.size();j++){
588 0 : if (handledEdges[j].index==edge){
589 : HandledEdge& hE=handledEdges[j];
590 : found=true;
591 : // set new from/to values
592 0 : for (size_t k=0;k<2;k++){
593 0 : if (hE.from){
594 0 : m_vSCVF[hE.scvfIndex+k].To=scvID[k];
595 : } else {
596 0 : m_vSCVF[hE.scvfIndex+k].From=scvID[k];
597 : }
598 : }
599 : // mark edge so associated scvf is not overwritten
600 0 : faceEdge[i]=deleted;
601 : break;
602 : }
603 : }
604 0 : if (found==true) continue;
605 : }
606 : HandledEdge hEdge;
607 0 : hEdge.index=edge;
608 0 : hEdge.scvfIndex=m_numSCVF;
609 : MathVector<worldDim> edgeMidGlob;
610 : MathVector<dim> edgeMidLoc;
611 0 : for (int d=0;d<worldDim;d++){
612 0 : edgeMidGlob[d] = 0.5 * (vCornerCoords[edgeCo[0]][d] + vCornerCoords[edgeCo[1]][d]);
613 0 : edgeMidLoc[d] = 0.5 * (m_rRefElem.corner(edgeCo[0])[d] + m_rRefElem.corner(edgeCo[1])[d]);
614 : }
615 0 : for (size_t j=0;j<2;j++){
616 0 : hEdge.associatedSCV[j] = scvID[j];
617 0 : if (from==face){
618 0 : m_vSCVF[m_numSCVF].From = scvID[j];
619 0 : m_vSCVF[m_numSCVF].To = to;
620 0 : hEdge.from=true;
621 : } else {
622 0 : m_vSCVF[m_numSCVF].From = from;
623 0 : m_vSCVF[m_numSCVF].To = scvID[j];
624 0 : hEdge.from=false;
625 : }
626 0 : m_vSCVF[m_numSCVF].Normal = normal;
627 0 : m_vSCVF[m_numSCVF].vGloPos[0] = vCornerCoords[edgeCo[j]];
628 : m_vSCVF[m_numSCVF].vLocPos[0] = m_rRefElem.corner(edgeCo[j]);
629 : m_vSCVF[m_numSCVF].vGloPos[1] = edgeMidGlob;
630 : m_vSCVF[m_numSCVF].vLocPos[1] = edgeMidLoc;
631 : m_vSCVF[m_numSCVF].vGloPos[2] = globalBary;
632 : m_vSCVF[m_numSCVF].vLocPos[2] = localBary;
633 0 : m_vSCVF[m_numSCVF].numSH = rTrialSpace.num_sh();
634 0 : AveragePositions(m_vSCVF[m_numSCVF].localIP, m_vSCVF[m_numSCVF].vLocPos, m_vSCVF[m_numSCVF].numCo);
635 0 : AveragePositions(m_vSCVF[m_numSCVF].globalIP, m_vSCVF[m_numSCVF].vGloPos, m_vSCVF[m_numSCVF].numCo);
636 0 : rTrialSpace.shapes(&(m_vSCVF[m_numSCVF].vShape[0]), m_vSCVF[m_numSCVF].local_ip());
637 0 : rTrialSpace.grads(&(m_vSCVF[m_numSCVF].vLocalGrad[0]), m_vSCVF[m_numSCVF].local_ip());
638 0 : m_numSCVF++;
639 : }
640 0 : handledEdges.push_back(hEdge);
641 : }
642 : // scvfs inside the face
643 : // insert remaining inner scvfs into positions of edge associated scvs
644 0 : for (size_t j=0;j<numFaceCo;j++){
645 : // replaces edge associated scvf
646 0 : size_t ii = faceEdge[j];
647 0 : if (ii==deleted){
648 0 : ii = m_numSCVF;
649 0 : m_numSCVF++;
650 : }
651 0 : if (numFaceCo==3){
652 : // compute inner scvfs in triangular case
653 0 : m_vSCVF[ii].From = m_numDofs+3;
654 0 : m_vSCVF[ii].To = m_numDofs+j;
655 : size_t ind = face;
656 0 : if (j>0) ind = m_numSCV+j-1;
657 : m_vSCVF[ii].vLocPos[0] = m_vSCV[ind].vLocPos[1];
658 : m_vSCVF[ii].vLocPos[1] = m_vSCV[ind].vLocPos[2];
659 : m_vSCVF[ii].vGloPos[0] = m_vSCV[ind].vGloPos[1];
660 : m_vSCVF[ii].vGloPos[1] = m_vSCV[ind].vGloPos[2];
661 : }else{
662 : // compute inner scvfs in quadrilateral case
663 0 : m_vSCVF[ii].To = m_numDofs+j;
664 0 : m_vSCVF[ii].From = m_numDofs + ((j+1) % 4);
665 0 : for (int d=0;d<worldDim;d++){
666 0 : m_vSCVF[ii].vLocPos[0][d] = 0.5*( m_rRefElem.corner(faceCo[j])[d] + m_rRefElem.corner(faceCo[(j+1) % 4])[d] );
667 0 : m_vSCVF[ii].vGloPos[0][d] = 0.5*( vCornerCoords[faceCo[j]][d] + vCornerCoords[faceCo[(j+1) % 4]][d] );
668 : }
669 : m_vSCVF[ii].vLocPos[1] = localMidpoint;
670 : m_vSCVF[ii].vGloPos[1] = globalMidpoint;
671 : }
672 : m_vSCVF[ii].vLocPos[2] = localBary;
673 : m_vSCVF[ii].vGloPos[2] = globalBary;
674 0 : m_vSCVF[ii].numSH = rTrialSpace.num_sh();
675 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[ii].Normal,m_vSCVF[ii].vGloPos);
676 0 : AveragePositions(m_vSCVF[ii].globalIP,m_vSCVF[ii].vGloPos,3);
677 0 : AveragePositions(m_vSCVF[ii].localIP,m_vSCVF[ii].vLocPos,3);
678 0 : rTrialSpace.shapes(&(m_vSCVF[ii].vShape[0]), m_vSCVF[ii].local_ip());
679 0 : rTrialSpace.grads(&(m_vSCVF[ii].vLocalGrad[0]), m_vSCVF[ii].local_ip());
680 : }
681 : // insert new constrained dof object
682 0 : m_vCD[m_numConstrainedDofs].i = face;
683 0 : m_vCD[m_numConstrainedDofs].numConstrainingDofs = 4;
684 0 : for (size_t i=0;i<4;i++){
685 0 : m_vCD[m_numConstrainedDofs].cDofInd[i] = m_numDofs+i;
686 : size_t ind=face;
687 0 : if (i>0) ind = m_numSCV+i-1;
688 0 : m_vCD[m_numConstrainedDofs].cDofWeights[i] = (number)m_vSCV[ind].Vol / faceVol;
689 : }
690 0 : m_numSCV+=3+keepOffset;
691 0 : m_numDofs+=4;
692 0 : m_numConstrainedDofs+=1;
693 0 : m_roid = ROID_UNKNOWN;
694 : }
695 0 : }// end of hanging node check
696 :
697 : //\todo compute with one virt. call
698 : // compute jacobian for linear mapping
699 0 : if(rMapping.is_linear())
700 : {
701 : MathMatrix<worldDim,dim> JtInv;
702 0 : rMapping.jacobian_transposed_inverse(JtInv, m_vSCVF[0].local_ip());
703 0 : const number detJ = rMapping.sqrt_gram_det(m_vSCVF[0].local_ip());
704 :
705 0 : for(size_t i = 0; i < num_scvf(); ++i)
706 : {
707 : m_vSCVF[i].JtInv = JtInv;
708 0 : m_vSCVF[i].detj = detJ;
709 : }
710 :
711 0 : for(size_t i = 0; i < num_scv(); ++i)
712 : {
713 : m_vSCV[i].JtInv = JtInv;
714 0 : m_vSCV[i].detj = detJ;
715 : }
716 : }
717 : // else compute jacobian for each integration point
718 : else
719 : {
720 0 : for(size_t i = 0; i < num_scvf(); ++i)
721 : {
722 0 : rMapping.jacobian_transposed_inverse(m_vSCVF[i].JtInv, m_vSCVF[i].local_ip());
723 0 : m_vSCVF[i].detj = rMapping.sqrt_gram_det(m_vSCVF[i].local_ip());
724 : }
725 0 : for(size_t i = 0; i < num_scv(); ++i)
726 : {
727 0 : rMapping.jacobian_transposed_inverse(m_vSCV[i].JtInv, m_vSCV[i].local_ip());
728 0 : m_vSCV[i].detj = rMapping.sqrt_gram_det(m_vSCV[i].local_ip());
729 : }
730 : }
731 :
732 : // compute global gradients
733 0 : for(size_t i = 0; i < num_scvf(); ++i)
734 0 : for(size_t sh = 0; sh < scvf(i).num_sh(); ++sh)
735 : MatVecMult(m_vSCVF[i].vGlobalGrad[sh], m_vSCVF[i].JtInv, m_vSCVF[i].vLocalGrad[sh]);
736 :
737 0 : for(size_t i = 0; i < num_scv(); ++i)
738 0 : for(size_t sh = 0; sh < scv(i).num_sh(); ++sh)
739 : MatVecMult(m_vSCV[i].vGlobalGrad[sh], m_vSCV[i].JtInv, m_vSCV[i].vLocalGrad[sh]);
740 :
741 : // copy ip points in list (SCVF)
742 0 : for(size_t i = 0; i < num_scvf(); ++i)
743 : m_vGlobSCVF_IP[i] = scvf(i).global_ip();
744 :
745 : // if no boundary subsets required, return
746 0 : if(num_boundary_subsets() == 0 || ish == NULL) return;
747 0 : else update_boundary_faces(pElem, vCornerCoords, ish);
748 : }
749 :
750 :
751 : /// update geometric data for given element (no shapes)
752 : template <int TDim, int TWorldDim>
753 0 : void DimCRFVGeometry<TDim, TWorldDim>::
754 : update_geometric_data(GridObject* pElem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
755 : {
756 : // If already update for this element, do nothing
757 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
758 :
759 : // refresh local data, if different roid given
760 0 : if(m_roid != pElem->reference_object_id())
761 : {
762 : // remember new roid
763 0 : m_roid = (ReferenceObjectID) pElem->reference_object_id();
764 :
765 : // update local data
766 0 : update_local_data();
767 : }
768 :
769 : // get reference element
770 : try{
771 : const DimReferenceElement<dim>& rRefElem
772 0 : = ReferenceElementProvider::get<dim>(m_roid);
773 :
774 : // compute barycenter coordinates
775 : globalBary = vCornerCoords[0];
776 0 : for (size_t j=1;j<rRefElem.num(0);j++){
777 0 : globalBary+=vCornerCoords[j];
778 : }
779 0 : globalBary*=1./(number)rRefElem.num(0);
780 :
781 : // compute global informations for scvf
782 0 : for(size_t i = 0; i < num_scvf(); ++i)
783 : {
784 0 : for (size_t j=0;j<m_vSCVF[i].numCo-1;j++){
785 0 : m_vSCVF[i].vGloPos[j]=vCornerCoords[rRefElem.id(dim-2,i,0,j)];
786 : }
787 : m_vSCVF[i].vGloPos[m_vSCVF[i].numCo-1]=globalBary;
788 0 : AveragePositions(m_vSCVF[i].globalIP, m_vSCVF[i].vGloPos, m_vSCVF[i].numCo);
789 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[i].Normal,m_vSCVF[i].vGloPos);// face_type0 identical to scvf type
790 : }
791 :
792 : // compute size of scv
793 0 : for(size_t i = 0; i < num_scv(); ++i)
794 : {
795 : // side nodes in reverse order to fulfill standard element order
796 0 : for (int j=0;j<m_vSCV[i].numCorners-1;j++){
797 0 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-2-j]=vCornerCoords[rRefElem.id(dim-1,i,0,j)];
798 : }
799 0 : AveragePositions(m_vGlobUnkCoords[i], m_vSCV[i].vGloPos, m_vSCV[i].numCorners-1);
800 : m_vSCV[i].vGlobIP = m_vGlobUnkCoords[i];
801 :
802 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-1]=globalBary;
803 : // compute volume of scv and normal to associated element face
804 : //CRSCVSizeAndNormal<dim>(m_vSCV[i].Vol,m_vSCV[i].Normal,m_vSCV[i].vGloPos,m_vSCV[i].numCorners);
805 0 : if (m_vSCV[i].numCorners-1==dim){
806 0 : m_vSCV[i].Vol = ElementSize<scv_type0,worldDim>(m_vSCV[i].vGloPos);
807 0 : ElementNormal<face_type0, worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
808 : } else { // m_vSCV[i].numCorners-2==dim , only possible in 3d (pyramid)
809 0 : m_vSCV[i].Vol = ElementSize<scv_type1,worldDim>(m_vSCV[i].vGloPos);
810 0 : ElementNormal<face_type1,worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
811 : };
812 : // nodes are in reverse order therefore reverse sign to get outward normal
813 : m_vSCV[i].Normal*=-1;
814 : }
815 :
816 : }
817 0 : UG_CATCH_THROW("DimCRFVGeometry: geometric update failed.");
818 : }
819 :
820 :
821 : template <int TDim, int TWorldDim>
822 0 : void DimCRFVGeometry<TDim, TWorldDim>::
823 : update_boundary_faces(GridObject* pElem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
824 : {
825 : // get grid
826 0 : Grid& grid = *(ish->grid());
827 :
828 : // vector of subset indices of side
829 : std::vector<int> vSubsetIndex;
830 :
831 : // get subset indices for sides (i.e. edge in 2d, faces in 3d)
832 : if(dim == 1) {
833 : std::vector<Vertex*> vVertex;
834 : CollectVertices(vVertex, grid, pElem);
835 : vSubsetIndex.resize(vVertex.size());
836 : for(size_t i = 0; i < vVertex.size(); ++i)
837 : vSubsetIndex[i] = ish->get_subset_index(vVertex[i]);
838 : }
839 : if(dim == 2) {
840 : std::vector<Edge*> vEdges;
841 0 : CollectEdgesSorted(vEdges, grid, pElem);
842 0 : vSubsetIndex.resize(vEdges.size());
843 0 : for(size_t i = 0; i < vEdges.size(); ++i)
844 0 : vSubsetIndex[i] = ish->get_subset_index(vEdges[i]);
845 0 : }
846 : if(dim == 3) {
847 : std::vector<Face*> vFaces;
848 0 : CollectFacesSorted(vFaces, grid, pElem);
849 0 : vSubsetIndex.resize(vFaces.size());
850 0 : for(size_t i = 0; i < vFaces.size(); ++i)
851 0 : vSubsetIndex[i] = ish->get_subset_index(vFaces[i]);
852 0 : }
853 :
854 : try{
855 : // const DimReferenceElement<dim>& rRefElem
856 : // = ReferenceElementProvider::get<dim>(m_roid);
857 :
858 0 : DimReferenceMapping<dim, worldDim>& rMapping = ReferenceMappingProvider::get<dim, worldDim>(m_roid);
859 0 : rMapping.update(vCornerCoords);
860 :
861 : const LocalShapeFunctionSet<dim>& TrialSpace =
862 0 : LocalFiniteElementProvider::get<dim>(m_roid, LFEID(LFEID::CROUZEIX_RAVIART, dim, 1));
863 :
864 : // loop requested subset
865 : typename std::map<int, std::vector<BF> >::iterator it;
866 0 : for (it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); ++it)
867 : {
868 : // get subset index
869 0 : const int bndIndex = (*it).first;
870 :
871 : // get vector of BF for element
872 0 : std::vector<BF>& vBF = (*it).second;
873 :
874 : // clear vector
875 : vBF.clear();
876 :
877 : // current number of bf
878 : size_t curr_bf = 0;
879 :
880 : // loop sides of element
881 0 : for(size_t side = 0; side < vSubsetIndex.size(); ++side)
882 : {
883 : // skip non boundary sides
884 0 : if(vSubsetIndex[side] != bndIndex) continue;
885 :
886 : // number of corners of side
887 : // const int coOfSide = rRefElem.num(dim-1, side, 0); todo use somewhere?
888 :
889 : // resize vector
890 0 : vBF.resize(curr_bf + 1);
891 :
892 : // fill BF with data from associated SCV
893 : BF& bf = vBF[curr_bf];
894 0 : bf.nodeID = m_vSCV[side].nodeID;
895 :
896 : bf.localIP = m_vSCV[side].vLocIP;
897 : bf.globalIP = m_vSCV[side].vGlobIP;
898 :
899 : bf.Normal = m_vSCV[side].Normal;
900 : // compute volume
901 0 : bf.Vol = VecTwoNorm(bf.Normal);
902 :
903 0 : bf.numCo = m_vSCV[side].numCorners-1;
904 :
905 : // compute shapes and grads
906 0 : bf.numSH = TrialSpace.num_sh();
907 0 : TrialSpace.shapes(&(bf.vShape[0]), bf.localIP);
908 0 : TrialSpace.grads(&(bf.vLocalGrad[0]), bf.localIP);
909 :
910 : // get reference mapping
911 0 : rMapping.jacobian_transposed_inverse(bf.JtInv, bf.localIP);
912 0 : bf.detj = rMapping.sqrt_gram_det(bf.localIP);
913 :
914 : // compute global gradients
915 0 : for(size_t sh = 0 ; sh < bf.num_sh(); ++sh)
916 : MatVecMult(bf.vGlobalGrad[sh], bf.JtInv, bf.vLocalGrad[sh]);
917 :
918 : // increase curr_bf
919 : ++curr_bf;
920 : }
921 : }
922 :
923 : }
924 0 : UG_CATCH_THROW("DimCRFVGeometry: update failed.");
925 0 : }
926 :
927 : ////////////////////////////////////////////////////////////////////////////////
928 : ////////////////////////////////////////////////////////////////////////////////
929 : //// Methods for CRFVGeometry class
930 : ////////////////////////////////////////////////////////////////////////////////
931 : ////////////////////////////////////////////////////////////////////////////////
932 :
933 : template < typename TElem, int TWorldDim>
934 0 : CRFVGeometry<TElem, TWorldDim>::
935 : CRFVGeometry()
936 0 : : m_pElem(NULL), m_rRefElem(Provider<ref_elem_type>::get()),
937 0 : m_rTrialSpace(Provider<local_shape_fct_set_type>::get())
938 : {
939 0 : update_local_data();
940 0 : }
941 :
942 : template < typename TElem, int TWorldDim>
943 0 : void CRFVGeometry<TElem, TWorldDim>::
944 : update_local_data()
945 : {
946 : // compute barycenter coordinates
947 0 : localBary = m_rRefElem.corner(0);
948 0 : for (size_t j=1;j<m_rRefElem.num(0);j++){
949 : localBary+=m_rRefElem.corner(j);
950 : }
951 0 : localBary*=1./(number)m_rRefElem.num(0);
952 :
953 : // set up local informations for SubControlVolumeFaces (scvf)
954 : // each scvf is associated to one vertex (2d) / edge (3d) of the element
955 0 : for(size_t i = 0; i < numSCVF; ++i)
956 : {
957 : // this scvf separates the given edges/faces
958 0 : m_vSCVF[i].From = m_rRefElem.id(dim-2, i, dim-1, 0);// todo handle dim==1
959 0 : m_vSCVF[i].To = m_rRefElem.id(dim-2, i, dim-1, 1);
960 :
961 0 : for (size_t j=0;j<m_vSCVF[i].numCo-1;j++){
962 0 : m_vSCVF[i].vLocPos[j]=m_rRefElem.corner(m_rRefElem.id(dim-2,i,0,j));
963 : }
964 :
965 : m_vSCVF[i].vLocPos[m_vSCVF[i].numCo-1]=localBary;
966 :
967 0 : AveragePositions(m_vSCVF[i].localIP, m_vSCVF[i].vLocPos, m_vSCVF[i].numCo);
968 : }
969 :
970 : // set up local informations for SubControlVolumes (scv)
971 : // each scv is associated to one edge (2d) / face (3d) of the element
972 0 : for(size_t i = 0; i < numSCV; ++i)
973 : {
974 : // store associated node
975 0 : m_vSCV[i].nodeID = i;
976 :
977 0 : m_vSCV[i].numCorners = m_rRefElem.num(dim-1,i,0)+1;
978 0 : for (int j=0;j<m_vSCV[i].numCorners-1;j++){
979 0 : m_vSCV[i].vLocPos[m_vSCV[i].numCorners-2-j]=m_rRefElem.corner(m_rRefElem.id(dim-1,i,0,j));
980 : }
981 0 : AveragePositions(m_vLocUnkCoords[i], m_vSCV[i].vLocPos, m_vSCV[i].numCorners-1);
982 : m_vSCV[i].vLocIP=m_vLocUnkCoords[i];
983 :
984 : m_vSCV[i].vLocPos[m_vSCV[i].numCorners-1]=localBary;
985 : }
986 :
987 : /////////////////////////
988 : // Shapes and Derivatives
989 : /////////////////////////
990 :
991 : // compute Shapes and Derivatives
992 0 : for(size_t i = 0; i < numSCVF; ++i)
993 : {
994 0 : m_rTrialSpace.shapes(&(m_vSCVF[i].vShape[0]), m_vSCVF[i].local_ip());
995 0 : m_rTrialSpace.grads(&(m_vSCVF[i].vLocalGrad[0]), m_vSCVF[i].local_ip());
996 : }
997 :
998 0 : for(size_t i = 0; i < numSCV; ++i)
999 : {
1000 0 : m_rTrialSpace.shapes(&(m_vSCV[i].vShape[0]), m_vSCV[i].local_ip());
1001 0 : m_rTrialSpace.grads(&(m_vSCV[i].vLocalGrad[0]), m_vSCV[i].local_ip());
1002 : }
1003 :
1004 : // copy ip positions in a list for Sub Control Volumes Faces (SCVF)
1005 0 : for(size_t i = 0; i < num_scvf(); ++i)
1006 : m_vLocSCVF_IP[i] = scvf(i).local_ip();
1007 0 : }
1008 :
1009 :
1010 : /// update data for given element
1011 : template < typename TElem, int TWorldDim>
1012 0 : void CRFVGeometry<TElem, TWorldDim>::
1013 : update(GridObject* elem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
1014 : {
1015 : UG_ASSERT(dynamic_cast<TElem*>(elem) != NULL, "Wrong element type.");
1016 : TElem* pElem = static_cast<TElem*>(elem);
1017 :
1018 : // If already update for this element, do nothing
1019 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
1020 :
1021 : // compute barycenter coordinates
1022 : globalBary = vCornerCoords[0];
1023 : m_vCo[0] = vCornerCoords[0];
1024 0 : for (size_t j=1;j<m_rRefElem.num(0);j++){
1025 0 : globalBary+=vCornerCoords[j];
1026 : m_vCo[j] = vCornerCoords[j];
1027 : }
1028 0 : globalBary*=1./(number)m_rRefElem.num(0);
1029 :
1030 : // compute global informations for scvf
1031 0 : for(size_t i = 0; i < num_scvf(); ++i)
1032 : {
1033 0 : for (size_t j=0;j<m_vSCVF[i].numCo-1;j++){
1034 0 : m_vSCVF[i].vGloPos[j]=vCornerCoords[m_rRefElem.id(dim-2,i,0,j)];
1035 : }
1036 : m_vSCVF[i].vGloPos[m_vSCVF[i].numCo-1]=globalBary;
1037 0 : AveragePositions(m_vSCVF[i].globalIP, m_vSCVF[i].vGloPos, m_vSCVF[i].numCo);
1038 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[i].Normal,m_vSCVF[i].vGloPos);// face_type0 identical to scvf type
1039 : }
1040 :
1041 : // compute size of scv
1042 0 : for(size_t i = 0; i < num_scv(); ++i)
1043 : {
1044 : // side nodes in reverse order to fulfill standard element order
1045 0 : for (int j=0;j<m_vSCV[i].numCorners-1;j++){
1046 0 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-2-j]=vCornerCoords[m_rRefElem.id(dim-1,i,0,j)];
1047 : }
1048 0 : AveragePositions(m_vGlobUnkCoords[i], m_vSCV[i].vGloPos, m_vSCV[i].numCorners-1);
1049 : m_vSCV[i].vGlobIP = m_vGlobUnkCoords[i];
1050 :
1051 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-1]=globalBary;
1052 : // compute volume of scv and normal to associated element face
1053 : //CRSCVSizeAndNormal<dim>(m_vSCV[i].Vol,m_vSCV[i].Normal,m_vSCV[i].vGloPos,m_vSCV[i].numCorners);
1054 0 : if (m_vSCV[i].numCorners-1==dim){
1055 0 : m_vSCV[i].Vol = ElementSize<scv_type0,worldDim>(m_vSCV[i].vGloPos);
1056 0 : ElementNormal<face_type0, worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
1057 : } else { // m_vSCV[i].numCorners-2==dim , only possible in 3d (pyramid)
1058 0 : m_vSCV[i].Vol = ElementSize<scv_type1,worldDim>(m_vSCV[i].vGloPos);
1059 0 : ElementNormal<face_type1, worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
1060 : };
1061 : // nodes are in reverse order therefore reverse sign to get outward normal
1062 : m_vSCV[i].Normal*=-1;
1063 : }
1064 :
1065 : // Shapes and Derivatives
1066 0 : m_mapping.update(vCornerCoords);
1067 :
1068 : // compute jacobian for linear mapping
1069 : if(ReferenceMapping<ref_elem_type, worldDim>::isLinear)
1070 : {
1071 : MathMatrix<worldDim,dim> JtInv;
1072 0 : m_mapping.jacobian_transposed_inverse(JtInv, m_vSCVF[0].local_ip());
1073 0 : const number detJ = m_mapping.sqrt_gram_det(m_vSCVF[0].local_ip());
1074 :
1075 0 : for(size_t i = 0; i < num_scvf(); ++i)
1076 : {
1077 : m_vSCVF[i].JtInv = JtInv;
1078 0 : m_vSCVF[i].detj = detJ;
1079 : }
1080 :
1081 0 : for(size_t i = 0; i < num_scv(); ++i)
1082 : {
1083 : m_vSCV[i].JtInv = JtInv;
1084 0 : m_vSCV[i].detj = detJ;
1085 : }
1086 : }
1087 : // else compute jacobian for each integration point
1088 : else
1089 : {
1090 0 : for(size_t i = 0; i < num_scvf(); ++i)
1091 : {
1092 0 : m_mapping.jacobian_transposed_inverse(m_vSCVF[i].JtInv, m_vSCVF[i].local_ip());
1093 0 : m_vSCVF[i].detj = m_mapping.sqrt_gram_det(m_vSCVF[i].local_ip());
1094 : }
1095 0 : for(size_t i = 0; i < num_scv(); ++i)
1096 : {
1097 0 : m_mapping.jacobian_transposed_inverse(m_vSCV[i].JtInv, m_vSCV[i].local_ip());
1098 0 : m_vSCV[i].detj = m_mapping.sqrt_gram_det(m_vSCV[i].local_ip());
1099 : }
1100 : }
1101 :
1102 : // compute global gradients
1103 0 : for(size_t i = 0; i < num_scvf(); ++i)
1104 0 : for(size_t sh = 0; sh < scvf(i).num_sh(); ++sh)
1105 : MatVecMult(m_vSCVF[i].vGlobalGrad[sh], m_vSCVF[i].JtInv, m_vSCVF[i].vLocalGrad[sh]);
1106 :
1107 0 : for(size_t i = 0; i < num_scv(); ++i)
1108 0 : for(size_t sh = 0; sh < scv(i).num_sh(); ++sh)
1109 : MatVecMult(m_vSCV[i].vGlobalGrad[sh], m_vSCV[i].JtInv, m_vSCV[i].vLocalGrad[sh]);
1110 :
1111 : // copy ip points in list (SCVF)
1112 0 : for(size_t i = 0; i < num_scvf(); ++i)
1113 : m_vGlobSCVF_IP[i] = scvf(i).global_ip();
1114 :
1115 : // if no boundary subsets required, return
1116 0 : if(num_boundary_subsets() == 0 || ish == NULL) return;
1117 0 : else update_boundary_faces(pElem, vCornerCoords, ish);
1118 : }
1119 :
1120 : template < typename TElem, int TWorldDim>
1121 0 : void CRFVGeometry<TElem, TWorldDim>::
1122 : update_boundary_faces(GridObject* pElem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
1123 : {
1124 : // get grid
1125 0 : Grid& grid = *(ish->grid());
1126 :
1127 : // vector of subset indices of side
1128 : std::vector<int> vSubsetIndex;
1129 :
1130 : // get subset indices for sides (i.e. edge in 2d, faces in 3d)
1131 : if(dim == 1) {
1132 : std::vector<Vertex*> vVertex;
1133 : CollectVertices(vVertex, grid, pElem);
1134 : vSubsetIndex.resize(vVertex.size());
1135 : for(size_t i = 0; i < vVertex.size(); ++i)
1136 : vSubsetIndex[i] = ish->get_subset_index(vVertex[i]);
1137 : }
1138 : if(dim == 2) {
1139 : std::vector<Edge*> vEdges;
1140 0 : CollectEdgesSorted(vEdges, grid, pElem);
1141 0 : vSubsetIndex.resize(vEdges.size());
1142 0 : for(size_t i = 0; i < vEdges.size(); ++i)
1143 0 : vSubsetIndex[i] = ish->get_subset_index(vEdges[i]);
1144 0 : }
1145 : if(dim == 3) {
1146 : std::vector<Face*> vFaces;
1147 0 : CollectFacesSorted(vFaces, grid, pElem);
1148 0 : vSubsetIndex.resize(vFaces.size());
1149 0 : for(size_t i = 0; i < vFaces.size(); ++i)
1150 0 : vSubsetIndex[i] = ish->get_subset_index(vFaces[i]);
1151 0 : }
1152 :
1153 : // loop requested subset
1154 : typename std::map<int, std::vector<BF> >::iterator it;
1155 0 : for (it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); ++it)
1156 : {
1157 : // get subset index
1158 0 : const int bndIndex = (*it).first;
1159 :
1160 : // get vector of BF for element
1161 0 : std::vector<BF>& vBF = (*it).second;
1162 :
1163 : // clear vector
1164 : vBF.clear();
1165 :
1166 : // current number of bf
1167 : size_t curr_bf = 0;
1168 :
1169 : // loop sides of element
1170 0 : for(size_t side = 0; side < vSubsetIndex.size(); ++side)
1171 : {
1172 : // skip non boundary sides
1173 0 : if(vSubsetIndex[side] != bndIndex) continue;
1174 :
1175 : // resize vector
1176 0 : vBF.resize(curr_bf + 1);
1177 :
1178 : // fill BF with data from associated SCV
1179 : BF& bf = vBF[curr_bf];
1180 0 : bf.nodeID = m_vSCV[side].nodeID;
1181 :
1182 : bf.localIP = m_vSCV[side].vLocIP;
1183 : bf.globalIP = m_vSCV[side].vGlobIP;
1184 :
1185 : bf.Normal = m_vSCV[side].Normal;
1186 : // compute volume
1187 0 : bf.Vol = VecTwoNorm(bf.Normal);
1188 :
1189 0 : bf.numCo = m_vSCV[side].numCorners-1;
1190 :
1191 : // compute shapes and grads
1192 0 : m_rTrialSpace.shapes(&(bf.vShape[0]), bf.localIP);
1193 0 : m_rTrialSpace.grads(&(bf.vLocalGrad[0]), bf.localIP);
1194 :
1195 : // get reference mapping
1196 0 : m_mapping.jacobian_transposed_inverse(bf.JtInv, bf.localIP);
1197 0 : bf.detj = m_mapping.sqrt_gram_det(bf.localIP);
1198 :
1199 : // compute global gradients
1200 0 : for(size_t sh = 0 ; sh < bf.num_sh(); ++sh)
1201 : MatVecMult(bf.vGlobalGrad[sh], bf.JtInv, bf.vLocalGrad[sh]);
1202 :
1203 : // increase curr_bf
1204 : ++curr_bf;
1205 : }
1206 : }
1207 0 : }
1208 :
1209 : #ifdef UG_DIM_2
1210 : template class DimCRFVGeometry<2, 2>;
1211 : template class CRFVGeometry<Triangle, 2>;
1212 : template class CRFVGeometry<Quadrilateral, 2>;
1213 : #endif
1214 :
1215 : #ifdef UG_DIM_3
1216 : template class DimCRFVGeometry<3, 3>;
1217 : template class CRFVGeometry<Triangle, 3>;
1218 : template class CRFVGeometry<Quadrilateral, 3>;
1219 :
1220 : template class CRFVGeometry<Tetrahedron, 3>;
1221 : template class CRFVGeometry<Prism, 3>;
1222 : template class CRFVGeometry<Pyramid, 3>;
1223 : template class CRFVGeometry<Hexahedron, 3>;
1224 : #endif
1225 :
1226 : } // end namespace ug
|
