Line data Source code
1 : /*
2 : * Copyright (c) 2013-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 "hfvcr_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 : //// Methods for HCRFVGeometry class
49 : ////////////////////////////////////////////////////////////////////////////////
50 : ////////////////////////////////////////////////////////////////////////////////
51 :
52 : template < typename TElem, int TWorldDim>
53 0 : HCRFVGeometry<TElem, TWorldDim>::
54 : HCRFVGeometry()
55 0 : : m_pElem(NULL), m_rRefElem(Provider<ref_elem_type>::get()),
56 0 : m_rTrialSpace(Provider<local_shape_fct_set_type>::get())
57 : {
58 0 : update_local_data();
59 0 : }
60 :
61 : template < typename TElem, int TWorldDim>
62 0 : void HCRFVGeometry<TElem, TWorldDim>::
63 : update_local_data()
64 : {
65 : // compute barycenter coordinates
66 0 : 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 edge of the element
74 0 : for(size_t i = 0; i < numNaturalSCVF; ++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 :
86 0 : AveragePositions(m_vSCVF[i].localIP, m_vSCVF[i].vLocPos, m_vSCVF[i].numCo);
87 : }
88 :
89 : // set up local informations for SubControlVolumes (scv)
90 : // each scv is associated to one corner of the element
91 0 : for(size_t i = 0; i < numNaturalSCV; ++i)
92 : {
93 : // store associated node
94 0 : m_vSCV[i].nodeID = i;
95 :
96 0 : m_vSCV[i].numCorners = m_rRefElem.num(dim-1,i,0)+1;
97 0 : for (int j=0;j<m_vSCV[i].numCorners-1;j++){
98 0 : m_vSCV[i].vLocPos[m_vSCV[i].numCorners-2-j]=m_rRefElem.corner(m_rRefElem.id(dim-1,i,0,j));
99 : }
100 0 : AveragePositions(m_vLocUnkCoords[i], m_vSCV[i].vLocPos, m_vSCV[i].numCorners-1);
101 : m_vSCV[i].vLocIP=m_vLocUnkCoords[i];
102 :
103 : m_vSCV[i].vLocPos[m_vSCV[i].numCorners-1]=localBary;
104 : }
105 :
106 : /////////////////////////
107 : // Shapes and Derivatives
108 : /////////////////////////
109 :
110 : // compute Shapes and Derivatives
111 0 : for(size_t i = 0; i < numNaturalSCVF; ++i)
112 : {
113 0 : m_rTrialSpace.shapes(&(m_vSCVF[i].vShape[0]), m_vSCVF[i].local_ip());
114 0 : m_rTrialSpace.grads(&(m_vSCVF[i].vLocalGrad[0]), m_vSCVF[i].local_ip());
115 : }
116 :
117 0 : for(size_t i = 0; i < numNaturalSCV; ++i)
118 : {
119 0 : m_rTrialSpace.shapes(&(m_vSCV[i].vShape[0]), m_vSCV[i].local_ip());
120 0 : m_rTrialSpace.grads(&(m_vSCV[i].vLocalGrad[0]), m_vSCV[i].local_ip());
121 : }
122 :
123 : // copy ip positions in a list for Sub Control Volumes Faces (SCVF)
124 0 : for(size_t i = 0; i < numNaturalSCVF; ++i)
125 : m_vLocSCVF_IP[i] = scvf(i).local_ip();
126 :
127 0 : localUpdateNecessary=false;
128 0 : numConstrainedDofs = 0;
129 0 : numSCV = numNaturalSCV;
130 0 : numSCVF = numNaturalSCVF;
131 0 : numDofs = numSCV;
132 0 : }
133 :
134 :
135 : /// update data for given element
136 : template < typename TElem, int TWorldDim>
137 0 : void HCRFVGeometry<TElem, TWorldDim>::
138 : update(GridObject* elem, const MathVector<worldDim>* vCornerCoords, const ISubsetHandler* ish)
139 : {
140 : UG_ASSERT(dynamic_cast<TElem*>(elem) != NULL, "Wrong element type.");
141 : TElem* pElem = static_cast<TElem*>(elem);
142 :
143 : // if already update for this element, do nothing
144 0 : if(m_pElem == pElem) return; else m_pElem = pElem;
145 :
146 : // get grid
147 0 : Grid& grid = *(ish->grid());
148 :
149 : // update local data if some of it has been overwritten by constrained object scv/scvf
150 0 : if (localUpdateNecessary) update_local_data();
151 :
152 : // compute barycenter coordinates
153 : globalBary = vCornerCoords[0];
154 : m_vCo[0] = vCornerCoords[0];
155 0 : for (size_t j=1;j<m_rRefElem.num(0);j++){
156 0 : globalBary+=vCornerCoords[j];
157 : m_vCo[j] = vCornerCoords[j];
158 : }
159 0 : globalBary*=1./(number)m_rRefElem.num(0);
160 :
161 : // compute global informations for scvf
162 0 : for(size_t i = 0; i < num_scvf(); ++i)
163 : {
164 0 : for (size_t j=0;j<m_vSCVF[i].numCo-1;j++){
165 0 : m_vSCVF[i].vGloPos[j]=vCornerCoords[m_rRefElem.id(dim-2,i,0,j)];
166 : }
167 : m_vSCVF[i].vGloPos[m_vSCVF[i].numCo-1]=globalBary;
168 0 : AveragePositions(m_vSCVF[i].globalIP, m_vSCVF[i].vGloPos, m_vSCVF[i].numCo);
169 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[i].Normal,m_vSCVF[i].vGloPos);// face_type0 identical to scvf type
170 : }
171 :
172 : // compute size of scv
173 0 : for(size_t i = 0; i < num_scv(); ++i)
174 : {
175 : // side nodes in reverse order to fulfill standard element order
176 0 : for (int j=0;j<m_vSCV[i].numCorners-1;j++){
177 0 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-2-j]=vCornerCoords[m_rRefElem.id(dim-1,i,0,j)];
178 : }
179 0 : AveragePositions(m_vGlobUnkCoords[i], m_vSCV[i].vGloPos, m_vSCV[i].numCorners-1);
180 : m_vSCV[i].vGlobIP = m_vGlobUnkCoords[i];
181 :
182 : m_vSCV[i].vGloPos[m_vSCV[i].numCorners-1]=globalBary;
183 : // compute volume of scv and normal to associated element face
184 : //CRSCVSizeAndNormal<dim>(m_vSCV[i].Vol,m_vSCV[i].Normal,m_vSCV[i].vGloPos,m_vSCV[i].numCorners);
185 0 : if (m_vSCV[i].numCorners-1==dim){
186 0 : m_vSCV[i].Vol = ElementSize<scv_type0,worldDim>(m_vSCV[i].vGloPos);
187 0 : ElementNormal<face_type0, worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
188 : } else { // m_vSCV[i].numCorners-2==dim , only possible in 3d (pyramid)
189 0 : m_vSCV[i].Vol = ElementSize<scv_type1,worldDim>(m_vSCV[i].vGloPos);
190 0 : ElementNormal<face_type1, worldDim>(m_vSCV[i].Normal, m_vSCV[i].vGloPos);
191 : };
192 : // nodes are in reverse order therefore reverse sign to get outward normal
193 : m_vSCV[i].Normal*=-1;
194 : }
195 :
196 : // check for hanging nodes
197 : if (dim==2){
198 : std::vector<Edge*> vEdges;
199 0 : CollectEdgesSorted(vEdges, grid, pElem);
200 0 : for(size_t side = 0; side < vEdges.size(); ++side){
201 0 : ConstrainingEdge* constrainingObj = dynamic_cast<ConstrainingEdge*>(vEdges[side]);
202 0 : if(constrainingObj == NULL) continue;
203 :
204 : // found constraining edge
205 : MathVector<worldDim> globalMidpoint = m_vSCV[side].vGlobIP;
206 : MathVector<dim> localMidpoint = m_vSCV[side].vLocIP;
207 : // get edge corners
208 : size_t edgeCo[2];
209 0 : for (size_t j=0;j<2;j++) edgeCo[j] = m_rRefElem.id(1,side,0,j);
210 : // compute dof positions on constraining edge,
211 : // replace dof "side" with first and insert second at the end
212 0 : for (int d=0;d<worldDim;d++)
213 0 : m_vGlobUnkCoords[side][d] = 0.5 * (vCornerCoords[edgeCo[0]][d] + globalMidpoint[d]);
214 0 : for (int d=0;d<dim;d++)
215 0 : m_vLocUnkCoords[side][d] = 0.5 * (m_rRefElem.corner(edgeCo[0])[d] + localMidpoint[d]);
216 0 : for (int d=0;d<worldDim;d++)
217 0 : m_vGlobUnkCoords[numSCV][d] = 0.5 * (vCornerCoords[edgeCo[1]][d] + globalMidpoint[d]);
218 0 : for (int d=0;d<dim;d++)
219 0 : m_vLocUnkCoords[numSCV][d] = 0.5 * (m_rRefElem.corner(edgeCo[1])[d] + localMidpoint[d]);
220 : // handle corresponding scvfs
221 0 : for (size_t j=0;j<2;j++){
222 0 : if (m_vSCVF[edgeCo[j]].From==side) m_vSCVF[edgeCo[j]].From=numDofs+j;
223 0 : if (m_vSCVF[edgeCo[j]].To==side) m_vSCVF[edgeCo[j]].To=numDofs+j;
224 : }
225 : // set up new scvs
226 : // scv "side" gets replaced by new one
227 : m_vSCV[side].Normal *= 0.5;
228 0 : m_vSCV[side].Vol *= 0.5;
229 0 : m_vSCV[side].nodeID = numDofs;
230 : m_vSCV[side].vGlobIP = m_vGlobUnkCoords[side];
231 : m_vSCV[side].vLocIP = m_vLocUnkCoords[side];
232 0 : m_rTrialSpace.shapes(&(m_vSCV[side].vShape[0]), m_vSCV[side].local_ip());
233 0 : m_rTrialSpace.grads(&(m_vSCV[side].vLocalGrad[0]), m_vSCV[side].local_ip());
234 : // second scv inserted at the end
235 0 : m_vSCV[numSCV].Normal = m_vSCV[side].Normal;
236 0 : m_vSCV[numSCV].Vol = m_vSCV[side].Vol;
237 0 : m_vSCV[numSCV].nodeID = numDofs+1;
238 : m_vSCV[numSCV].vGlobIP = m_vGlobUnkCoords[numSCV];
239 : m_vSCV[numSCV].vLocIP = m_vLocUnkCoords[numSCV];
240 0 : m_rTrialSpace.shapes(&(m_vSCV[numSCV].vShape[0]), m_vSCV[numSCV].local_ip());
241 0 : m_rTrialSpace.grads(&(m_vSCV[numSCV].vLocalGrad[0]), m_vSCV[numSCV].local_ip());
242 : // insert new scvf
243 0 : m_vSCVF[numSCVF].From = numDofs;
244 0 : m_vSCVF[numSCVF].To = numDofs+1;
245 : m_vSCVF[numSCVF].vLocPos[0] = localMidpoint;
246 : m_vSCVF[numSCVF].vLocPos[1] = localBary;
247 : m_vSCVF[numSCVF].vGloPos[0] = globalMidpoint;
248 : m_vSCVF[numSCVF].vGloPos[1] = globalBary;
249 0 : for (int d=0;d<dim;d++) m_vSCVF[numSCVF].localIP[d] = 0.5*(localMidpoint[d] + localBary[d]);
250 0 : for (int d=0;d<worldDim;d++) m_vSCVF[numSCVF].globalIP[d] = 0.5*(globalMidpoint[d] + globalBary[d]);
251 0 : m_rTrialSpace.shapes(&(m_vSCVF[numSCVF].vShape[0]), m_vSCVF[numSCVF].local_ip());
252 0 : m_rTrialSpace.grads(&(m_vSCVF[numSCVF].vLocalGrad[0]), m_vSCVF[numSCVF].local_ip());
253 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[numSCVF].Normal,m_vSCVF[numSCVF].vGloPos);
254 : // insert new constrained dof object
255 0 : m_vCD[numConstrainedDofs].i = side;
256 0 : m_vCD[numConstrainedDofs].numConstrainingDofs = 2;
257 0 : m_vCD[numConstrainedDofs].cDofInd[0] = numDofs;
258 0 : m_vCD[numConstrainedDofs].cDofInd[1] = numDofs+1;
259 0 : m_vCD[numConstrainedDofs].cDofWeights[0] = 0.5;
260 0 : m_vCD[numConstrainedDofs].cDofWeights[1] = 0.5;
261 0 : numSCV+=1;
262 0 : numSCVF+=1;
263 0 : numDofs+=2;
264 0 : numConstrainedDofs+=1;
265 0 : localUpdateNecessary = true;
266 : }
267 0 : } else {
268 : // dim = 3
269 : std::vector<Face*> vFaces;
270 0 : CollectFacesSorted(vFaces, grid, pElem);
271 : handledEdges.clear();
272 0 : for(size_t face = 0; face < vFaces.size(); ++face){
273 0 : ConstrainingFace* constrainingObj = dynamic_cast<ConstrainingFace*>(vFaces[face]);
274 0 : if(constrainingObj == NULL) continue;
275 : // found constraining face
276 : MathVector<worldDim> globalMidpoint = m_vSCV[face].vGlobIP;
277 : MathVector<dim> localMidpoint = m_vSCV[face].vLocIP;
278 0 : number faceVol = m_vSCV[face].Vol;
279 : MathVector<worldDim> faceNormal = m_vSCV[face].Normal;
280 : // get face corners and edges
281 : size_t faceCo[4];
282 : size_t faceEdge[4];
283 : // number of corners of face = number of edges
284 0 : size_t numFaceCo = m_rRefElem.num(2,face,0);
285 0 : for (size_t j=0;j<numFaceCo;j++) faceCo[j] = m_rRefElem.id(2,face,0,j);
286 0 : for (size_t j=0;j<numFaceCo;j++) faceEdge[j] = m_rRefElem.id(2,face,1,j);
287 : // compute coordinates of each face and fill scv values
288 0 : for (size_t i=0;i<numFaceCo;i++){
289 0 : size_t co = faceCo[i];
290 : size_t nOfEdges=0;
291 : size_t nbEdges[2];
292 : // find 2 edges in face belonging to node
293 0 : for (size_t j=0;j<m_rRefElem.num(0,co,1);j++){
294 0 : size_t candidate = m_rRefElem.id(0,co,1,j);
295 : bool found = false;
296 0 : for (size_t k=0;k<numFaceCo;k++){
297 0 : if (faceEdge[k]==candidate){
298 : found = true;
299 : break;
300 : }
301 : }
302 0 : if (found==true){
303 0 : nbEdges[nOfEdges] = candidate;
304 0 : nOfEdges++;
305 0 : if (nOfEdges==2) break;
306 : }
307 : }
308 : // in triangular case switch edges if necessary for correct orientation
309 0 : if (numFaceCo==3){
310 0 : if (faceEdge[i]==nbEdges[1]){
311 0 : nbEdges[1] = nbEdges[0];
312 0 : nbEdges[0] = faceEdge[i];
313 : }
314 : }
315 : // first new scv replaces scv nr face
316 : size_t ind = face;
317 : // others are inserted at the end
318 0 : if (i>0) ind = numSCV+i-1;
319 0 : m_vSCV[ind].vGloPos[0] = vCornerCoords[co];
320 : m_vSCV[ind].vLocPos[0] = m_rRefElem.corner(co);
321 0 : for (int d=0;d<worldDim;d++){
322 : // edge 0 midpoint
323 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] );
324 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] );
325 : // edge 1 midpoint
326 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] );
327 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] );
328 : }
329 0 : if (numFaceCo==4) m_vSCV[ind].vGloPos[2] = globalMidpoint;
330 : m_vSCV[ind].vGloPos[numFaceCo] = globalBary;
331 : m_vSCV[ind].vLocPos[numFaceCo] = localBary;
332 0 : m_vSCV[ind].numCorners = numFaceCo + 1;
333 0 : AveragePositions(m_vGlobUnkCoords[ind], m_vSCV[ind].vGloPos, m_vSCV[ind].numCorners-1);
334 0 : AveragePositions(m_vLocUnkCoords[ind], m_vSCV[ind].vLocPos, m_vSCV[ind].numCorners-1);
335 : m_vSCV[ind].vLocIP = m_vLocUnkCoords[ind];
336 : m_vSCV[ind].vGlobIP = m_vGlobUnkCoords[ind];
337 0 : if (numFaceCo==3) m_vSCV[ind].Vol = ElementSize<scv_type0,worldDim>(m_vSCV[ind].vGloPos);
338 0 : else m_vSCV[ind].Vol = ElementSize<scv_type1,worldDim>(m_vSCV[ind].vGloPos);
339 0 : if (m_vSCV[ind].Vol<0) m_vSCV[ind].Vol *= -1;
340 : m_vSCV[ind].Normal = faceNormal;
341 0 : m_vSCV[ind].Normal *= (number) m_vSCV[ind].Vol / faceVol;
342 0 : m_rTrialSpace.shapes(&(m_vSCV[ind].vShape[0]), m_vSCV[ind].local_ip());
343 0 : m_rTrialSpace.grads(&(m_vSCV[ind].vLocalGrad[0]), m_vSCV[ind].local_ip());
344 0 : m_vSCV[ind].nodeID = numDofs+i;
345 : }
346 : // compute inner scv in triangular case
347 0 : if (numFaceCo==3){
348 0 : number volSum=m_vSCV[face].Vol;
349 0 : for (size_t j=0;j<2;j++){
350 0 : volSum+=m_vSCV[numSCV+j].Vol;
351 : }
352 0 : size_t ind = numSCV+2;
353 0 : m_vSCV[ind].Vol = faceVol - volSum;
354 0 : m_vSCV[ind].nodeID=numDofs+3;
355 : m_vSCV[ind].Normal = faceNormal;
356 0 : m_vSCV[ind].Normal *= (number) m_vSCV[ind].Vol / faceVol;
357 : m_vSCV[ind].vGlobIP = m_vSCV[face].vGloPos[1];
358 : m_vSCV[ind].vLocIP = m_vSCV[face].vLocPos[1];
359 0 : for (size_t j=0;j<2;j++){
360 0 : m_vSCV[ind].vGlobIP += m_vSCV[numSCV+j].vGloPos[1];
361 : m_vSCV[ind].vLocIP += m_vSCV[numSCV+j].vLocPos[1];
362 : }
363 : m_vSCV[ind].vGlobIP *= (number)1.0/3.0;
364 : m_vSCV[ind].vLocIP *= (number)1.0/3.0;
365 : m_vGlobUnkCoords[ind] = m_vSCV[ind].vGlobIP;
366 : m_vLocUnkCoords[ind] = m_vSCV[ind].vLocIP;
367 0 : m_vSCV[ind].numCorners = numFaceCo + 1;
368 0 : m_rTrialSpace.shapes(&(m_vSCV[ind].vShape[0]), m_vSCV[ind].local_ip());
369 0 : m_rTrialSpace.grads(&(m_vSCV[ind].vLocalGrad[0]), m_vSCV[ind].local_ip());
370 : }
371 : // insert new scvfs, first the ones associated to edges of face
372 0 : for (size_t i=0;i<numFaceCo;i++){
373 0 : size_t edge = faceEdge[i];
374 0 : size_t from = m_vSCVF[edge].From;
375 0 : size_t to = m_vSCVF[edge].To;
376 : MathVector<worldDim> normal = m_vSCVF[edge].Normal;
377 : normal*=0.5;
378 : size_t edgeCo[2];
379 : size_t scvID[2];
380 0 : for (size_t j=0;j<2;j++){
381 0 : edgeCo[j] = m_rRefElem.id(1,edge,0,j);
382 : // find corresponding face vertex (= corresponding scv index)
383 0 : for (size_t k=0;k<numFaceCo;k++){
384 0 : if (faceCo[k]==edgeCo[j]){
385 0 : scvID[j] = numDofs+k;
386 0 : break;
387 : }
388 : }
389 : }
390 : // look if edge has already been handled
391 0 : if (numConstrainedDofs>0){
392 : bool found=false;
393 0 : for (size_t j=0;j<handledEdges.size();j++){
394 0 : if (handledEdges[j].index==edge){
395 : HandledEdge& hE=handledEdges[j];
396 : found=true;
397 : // set new from/to values
398 0 : for (size_t k=0;k<2;k++){
399 0 : if (hE.from){
400 0 : m_vSCVF[hE.scvfIndex+k].To=scvID[k];
401 : } else {
402 0 : m_vSCVF[hE.scvfIndex+k].From=scvID[k];
403 : }
404 : }
405 : // mark edge so associated scvf is not overwritten
406 0 : faceEdge[i]=deleted;
407 : break;
408 : }
409 : }
410 0 : if (found==true) continue;
411 : }
412 : HandledEdge hEdge;
413 0 : hEdge.index=edge;
414 0 : hEdge.scvfIndex=numSCVF;
415 : MathVector<worldDim> edgeMidGlob;
416 : MathVector<dim> edgeMidLoc;
417 0 : for (int d=0;d<worldDim;d++){
418 0 : edgeMidGlob[d] = 0.5 * (vCornerCoords[edgeCo[0]][d] + vCornerCoords[edgeCo[1]][d]);
419 0 : edgeMidLoc[d] = 0.5 * (m_rRefElem.corner(edgeCo[0])[d] + m_rRefElem.corner(edgeCo[1])[d]);
420 : }
421 0 : for (size_t j=0;j<2;j++){
422 0 : hEdge.associatedSCV[j] = scvID[j];
423 0 : if (from==face){
424 0 : m_vSCVF[numSCVF].From = scvID[j];
425 0 : m_vSCVF[numSCVF].To = to;
426 0 : hEdge.from=true;
427 : } else {
428 0 : m_vSCVF[numSCVF].From = from;
429 0 : m_vSCVF[numSCVF].To = scvID[j];
430 0 : hEdge.from=false;
431 : }
432 0 : m_vSCVF[numSCVF].Normal = normal;
433 0 : m_vSCVF[numSCVF].vGloPos[0] = vCornerCoords[edgeCo[j]];
434 0 : m_vSCVF[numSCVF].vLocPos[0] = m_rRefElem.corner(edgeCo[j]);
435 : m_vSCVF[numSCVF].vGloPos[1] = edgeMidGlob;
436 : m_vSCVF[numSCVF].vLocPos[1] = edgeMidLoc;
437 : m_vSCVF[numSCVF].vGloPos[2] = globalBary;
438 : m_vSCVF[numSCVF].vLocPos[2] = localBary;
439 0 : AveragePositions(m_vSCVF[numSCVF].localIP, m_vSCVF[numSCVF].vLocPos, m_vSCVF[numSCVF].numCo);
440 0 : AveragePositions(m_vSCVF[numSCVF].globalIP, m_vSCVF[numSCVF].vGloPos, m_vSCVF[numSCVF].numCo);
441 0 : m_rTrialSpace.shapes(&(m_vSCVF[numSCVF].vShape[0]), m_vSCVF[numSCVF].local_ip());
442 0 : m_rTrialSpace.grads(&(m_vSCVF[numSCVF].vLocalGrad[0]), m_vSCVF[numSCVF].local_ip());
443 0 : numSCVF++;
444 : }
445 0 : handledEdges.push_back(hEdge);
446 : }
447 : // scvfs inside the face
448 : // insert remaining inner scvfs into positions of edge associated scvs
449 0 : for (size_t j=0;j<numFaceCo;j++){
450 : // replaces edge associated scvf
451 0 : size_t ii = faceEdge[j];
452 0 : if (ii==deleted){
453 0 : ii = numSCVF;
454 0 : numSCVF++;
455 : }
456 0 : if (numFaceCo==3){
457 : // compute inner scvfs in triangular case
458 0 : m_vSCVF[ii].From = numDofs+3;
459 0 : m_vSCVF[ii].To = numDofs+j;
460 : size_t ind = face;
461 0 : if (j>0) ind = numSCV+j-1;
462 : m_vSCVF[ii].vLocPos[0] = m_vSCV[ind].vLocPos[1];
463 : m_vSCVF[ii].vLocPos[1] = m_vSCV[ind].vLocPos[2];
464 : m_vSCVF[ii].vGloPos[0] = m_vSCV[ind].vGloPos[1];
465 : m_vSCVF[ii].vGloPos[1] = m_vSCV[ind].vGloPos[2];
466 : }else{
467 : // compute inner scvfs in quadrilateral case
468 0 : m_vSCVF[ii].To = numDofs+j;
469 0 : m_vSCVF[ii].From = numDofs + ((j+1) % 4);
470 0 : for (int d=0;d<worldDim;d++){
471 0 : m_vSCVF[ii].vLocPos[0][d] = 0.5*( m_rRefElem.corner(faceCo[j])[d] + m_rRefElem.corner(faceCo[(j+1) % 4])[d] );
472 0 : m_vSCVF[ii].vGloPos[0][d] = 0.5*( vCornerCoords[faceCo[j]][d] + vCornerCoords[faceCo[(j+1) % 4]][d] );
473 : }
474 : m_vSCVF[ii].vLocPos[1] = localMidpoint;
475 : m_vSCVF[ii].vGloPos[1] = globalMidpoint;
476 : }
477 : m_vSCVF[ii].vLocPos[2] = localBary;
478 : m_vSCVF[ii].vGloPos[2] = globalBary;
479 0 : ElementNormal<face_type0,worldDim>(m_vSCVF[ii].Normal,m_vSCVF[ii].vGloPos);
480 0 : AveragePositions(m_vSCVF[ii].globalIP,m_vSCVF[ii].vGloPos,3);
481 0 : AveragePositions(m_vSCVF[ii].localIP,m_vSCVF[ii].vLocPos,3);
482 0 : m_rTrialSpace.shapes(&(m_vSCVF[ii].vShape[0]), m_vSCVF[ii].local_ip());
483 0 : m_rTrialSpace.grads(&(m_vSCVF[ii].vLocalGrad[0]), m_vSCVF[ii].local_ip());
484 : }
485 : // insert new constrained dof object
486 0 : m_vCD[numConstrainedDofs].i = face;
487 0 : m_vCD[numConstrainedDofs].numConstrainingDofs = 4;
488 0 : for (size_t i=0;i<4;i++){
489 0 : m_vCD[numConstrainedDofs].cDofInd[i] = numDofs+i;
490 : size_t ind=face;
491 0 : if (i>0) ind = numSCV+i-1;
492 0 : m_vCD[numConstrainedDofs].cDofWeights[i] = (number)m_vSCV[ind].Vol / faceVol;
493 : }
494 0 : numSCV+=3;
495 0 : numDofs+=4;
496 0 : numConstrainedDofs+=1;
497 0 : localUpdateNecessary = true;
498 : }
499 0 : }
500 :
501 :
502 : // Shapes and Derivatives
503 0 : m_mapping.update(vCornerCoords);
504 :
505 : // compute jacobian for linear mapping
506 : if(ReferenceMapping<ref_elem_type, worldDim>::isLinear)
507 : {
508 : MathMatrix<worldDim,dim> JtInv;
509 0 : m_mapping.jacobian_transposed_inverse(JtInv, m_vSCVF[0].local_ip());
510 0 : const number detJ = m_mapping.sqrt_gram_det(m_vSCVF[0].local_ip());
511 :
512 0 : for(size_t i = 0; i < num_scvf(); ++i)
513 : {
514 : m_vSCVF[i].JtInv = JtInv;
515 0 : m_vSCVF[i].detj = detJ;
516 : }
517 :
518 0 : for(size_t i = 0; i < num_scv(); ++i)
519 : {
520 : m_vSCV[i].JtInv = JtInv;
521 0 : m_vSCV[i].detj = detJ;
522 : }
523 : }
524 : // else compute jacobian for each integration point
525 : else
526 : {
527 0 : for(size_t i = 0; i < num_scvf(); ++i)
528 : {
529 0 : m_mapping.jacobian_transposed_inverse(m_vSCVF[i].JtInv, m_vSCVF[i].local_ip());
530 0 : m_vSCVF[i].detj = m_mapping.sqrt_gram_det(m_vSCVF[i].local_ip());
531 : }
532 0 : for(size_t i = 0; i < num_scv(); ++i)
533 : {
534 0 : m_mapping.jacobian_transposed_inverse(m_vSCV[i].JtInv, m_vSCV[i].local_ip());
535 0 : m_vSCV[i].detj = m_mapping.sqrt_gram_det(m_vSCV[i].local_ip());
536 : }
537 : }
538 :
539 : // compute global gradients
540 0 : for(size_t i = 0; i < num_scvf(); ++i)
541 0 : for(size_t sh = 0; sh < scvf(i).num_sh(); ++sh)
542 : MatVecMult(m_vSCVF[i].vGlobalGrad[sh], m_vSCVF[i].JtInv, m_vSCVF[i].vLocalGrad[sh]);
543 :
544 0 : for(size_t i = 0; i < num_scv(); ++i)
545 0 : for(size_t sh = 0; sh < scv(i).num_sh(); ++sh)
546 : MatVecMult(m_vSCV[i].vGlobalGrad[sh], m_vSCV[i].JtInv, m_vSCV[i].vLocalGrad[sh]);
547 :
548 : // copy ip points in list (SCVF)
549 0 : for(size_t i = 0; i < num_scvf(); ++i)
550 : m_vGlobSCVF_IP[i] = scvf(i).global_ip();
551 :
552 : // debug output
553 : // if (localUpdateNecessary==true) print();
554 :
555 : }
556 :
557 :
558 : // debug output
559 : template < typename TElem, int TWorldDim>
560 0 : void HCRFVGeometry<TElem, TWorldDim>::print()
561 : {
562 : UG_LOG("\nFVG hanging debug output\n");
563 0 : for(size_t i = 0; i < numSCV; ++i)
564 : {
565 : UG_LOG(i<<" SCV: ");
566 : UG_LOG("node_id=" << m_vSCV[i].node_id());
567 0 : UG_LOG(", local_pos="<< m_vSCV[i].local_ip());
568 0 : UG_LOG(", global_pos="<< m_vSCV[i].global_ip());
569 : UG_LOG(", vol=" << m_vSCV[i].volume());
570 : UG_LOG("\n");
571 : }
572 : UG_LOG("\n");
573 0 : for(size_t i = 0; i < numSCVF; ++i)
574 : {
575 : UG_LOG(i<<" SCVF: ");
576 : UG_LOG("from=" << m_vSCVF[i].from()<<", to="<<m_vSCVF[i].to());
577 0 : UG_LOG(", local_pos="<< m_vSCVF[i].local_ip());
578 0 : UG_LOG(", global_pos="<< m_vSCVF[i].global_ip());
579 0 : UG_LOG(", normal=" << m_vSCVF[i].normal());
580 : UG_LOG("\n Shapes:\n");
581 0 : for(size_t sh=0; sh < m_vSCVF[i].num_sh(); ++sh)
582 : {
583 : UG_LOG(" " <<sh << ": shape="<<m_vSCVF[i].shape(sh));
584 0 : UG_LOG(", global_grad="<<m_vSCVF[i].global_grad(sh));
585 0 : UG_LOG(", local_grad="<<m_vSCVF[i].local_grad(sh));
586 : UG_LOG("\n");
587 : }
588 : }
589 : UG_LOG("\n");
590 : UG_LOG("constrained dof indices:\n");
591 0 : for (size_t i=0;i<numConstrainedDofs;i++){
592 0 : UG_LOG("constrained index = " << m_vCD[i].i << "\n");
593 : UG_LOG("constraining indices: ");
594 0 : for (size_t j=0;j< m_vCD[i].numConstrainingDofs;j++) UG_LOG(m_vCD[i].cDofInd[j] << " ");
595 : UG_LOG("\n");
596 : UG_LOG("weights: ");
597 0 : for (size_t j=0;j< m_vCD[i].numConstrainingDofs;j++) UG_LOG(m_vCD[i].cDofWeights[j] << " ");
598 : UG_LOG("\n");
599 : }
600 : UG_LOG("\n");
601 0 : }
602 :
603 : #ifdef UG_DIM_2
604 : template class HCRFVGeometry<Triangle, 2>;
605 : template class HCRFVGeometry<Quadrilateral, 2>;
606 : #endif
607 :
608 : #ifdef UG_DIM_3
609 : template class HCRFVGeometry<Triangle, 3>;
610 : template class HCRFVGeometry<Quadrilateral, 3>;
611 : template class HCRFVGeometry<Tetrahedron, 3>;
612 : template class HCRFVGeometry<Prism, 3>;
613 : template class HCRFVGeometry<Pyramid, 3>;
614 : template class HCRFVGeometry<Hexahedron, 3>;
615 : #endif
616 :
617 : } // end namespace ug
|
