GCC Code Coverage Report

Directory: ./
File: FiniteElement/currentFiniteElementWithBd.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 94 131 71.8%
Branches: 100 196 51.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: JJ-F. Gerbeaun & V.Martin
13 //
14
15 // System includes
16
17 // External includes
18 #include <boost/numeric/ublas/matrix_proxy.hpp>
19
20 // Project includes
21 #include "currentFiniteElementWithBd.hpp"
22
23 namespace felisce {
24
25 16912 CurrentFiniteElementWithBd::CurrentFiniteElementWithBd(
26 const RefElement& refEle,
27 const GeoElement& geoEle,
28 const std::vector<DegreeOfExactness>& degOfExactness,
29 const DegreeOfExactness& degOfExactnessBd
30 16912 ): CurrentFiniteElement(refEle,geoEle,degOfExactness),
31 16912 m_numBdEle( geoEle.numBdEle() )
32 {
33
34
6/14
✓ Branch 1 taken 3402 times.
✓ Branch 2 taken 13510 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3402 times.
✓ Branch 7 taken 13510 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 13510 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 16912 times.
 16912 FEL_ASSERT( ( m_numBdEle == geoEle.numFace() && numRefCoor() == 3 ) || //3D type of element (tetra, prisms..)
35 ( m_numBdEle == geoEle.numEdge() && numRefCoor() == 2 ) || //2D (tria,quadrangles)
36 ( m_numBdEle == 2 && numRefCoor() == 1 ) ); //1D (segment)
37 //--------------------
38 /// 1/ build the curvilinear bd elements (for the measure, normal of the faces)
39 //--------------------
40
2/4
✓ Branch 0 taken 16912 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 16912 times.
✗ Branch 5 not taken.
 16912 m_bdEle = new CurvilinearFiniteElement* [m_numBdEle];
41
42 /// NOTE: I construct a CurvilinearFE for *each* bd element of the element
43 /// we could also define m_numBdEle=1 in the constructor, if only 1 face is needed... vm 07/2011
44
2/2
✓ Branch 0 taken 56258 times.
✓ Branch 1 taken 16912 times.
 73170 for(int ibd=0; ibd<m_numBdEle; ibd++ ) {
45
1/2
✓ Branch 2 taken 56258 times.
✗ Branch 3 not taken.
 112516 m_bdEle[ibd] = new CurvilinearFiniteElement( CurrentFiniteElement::refEle().boundaryRefElement(ibd),
46 56258 CurrentFiniteElement::geoEle().boundaryGeoElement(ibd),
47
2/4
✓ Branch 1 taken 56258 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 56258 times.
✗ Branch 5 not taken.
 56258 degOfExactnessBd );
48 }
49
50 //--------------------
51 /// 2/ compute the boundary quadrature points on the reference element (from the bd ref element)
52 //--------------------
53 /// coordinates of the points on the Bd element
54
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 UBlasMatrix pointCoorOnBd; //one could also save this info, if necessary...
55
56
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_pointBdEle.resize(m_numBdEle);
57
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_numPointBdEle.resize(m_numBdEle);
58
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_refQuadPointBdEle.resize(m_numBdEle);
59
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_numQuadPointBdEle.resize(m_numBdEle);
60
61 /// compute the image of the quadrature points (that live on the boundary ref element)
62 /// onto the boundary elements of the reference element.
63 /// (eg. : image of the quad points of a ref triangle onto the triangles of the ref tetra)
64
2/2
✓ Branch 0 taken 56258 times.
✓ Branch 1 taken 16912 times.
 73170 for ( int ibd = 0; ibd <m_numBdEle; ibd++ ) {
65
1/2
✓ Branch 1 taken 56258 times.
✗ Branch 2 not taken.
 56258 m_numPointBdEle[ibd] = bdEle(ibd).geoEle().numPoint();
66
67 // Get the points of the boundary element from the reference element
68
1/2
✓ Branch 3 taken 56258 times.
✗ Branch 4 not taken.
 56258 m_pointBdEle[ibd].resize( m_numPointBdEle[ibd] );
69
1/2
✓ Branch 2 taken 56258 times.
✗ Branch 3 not taken.
 56258 pointCoorOnBd.resize(m_numPointBdEle[ibd],m_numCoor);
70
2/2
✓ Branch 1 taken 138124 times.
✓ Branch 2 taken 56258 times.
 194382 for ( int iptBd = 0; iptBd < m_numPointBdEle[ibd] ; iptBd++ ) {
71
1/2
✓ Branch 2 taken 138124 times.
✗ Branch 3 not taken.
 138124 m_pointBdEle[ibd][ iptBd ] = CurrentFiniteElement::geoEle().pointOfBdEle( ibd, iptBd );
72
73 // update the coordinates of the reference element
74
2/2
✓ Branch 0 taken 332912 times.
✓ Branch 1 taken 138124 times.
 471036 for(int icoor=0; icoor<m_numCoor; icoor++) {
75
2/4
✓ Branch 4 taken 332912 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 332912 times.
✗ Branch 8 not taken.
 332912 pointCoorOnBd(iptBd,icoor) = CurrentFiniteElement::geoEle().pointCoor(m_pointBdEle[ibd][iptBd],icoor);
76 }
77 }
78
79 /*
80 std::cout << "pointBdEle (" << m_pointBdEle[ibd].size() << ")" << std::endl << "[";
81 for ( int iptBd = 0; iptBd < m_numPointBdEle[ibd] ;iptBd++ ) {std::cout << m_pointBdEle[ibd][iptBd] << ";";}
82 std::cout << "]" << std::endl;
83 std::cout << "pointCoorOnBd" << std::endl;
84 std::cout << pointCoorOnBd << std::endl;
85 */
86
87
1/2
✓ Branch 1 taken 56258 times.
✗ Branch 2 not taken.
 56258 m_numQuadPointBdEle[ibd] = bdEle(ibd).numQuadraturePoint();
88 //std::cout << "numQuadraturePointBd=" << m_numQuadPointBdEle[ibd] <<std::endl;
89
90
1/2
✓ Branch 3 taken 56258 times.
✗ Branch 4 not taken.
 56258 std::vector<UBlasVector> phiGeoBd(m_numQuadPointBdEle[ibd]);
91
1/2
✓ Branch 3 taken 56258 times.
✗ Branch 4 not taken.
 56258 m_refQuadPointBdEle[ibd] = std::vector<Point>(m_numQuadPointBdEle[ibd]);
92
93
2/2
✓ Branch 1 taken 156388 times.
✓ Branch 2 taken 56258 times.
 212646 for(int ig=0; ig< m_numQuadPointBdEle[ibd]; ig++) {
94
1/2
✓ Branch 3 taken 156388 times.
✗ Branch 4 not taken.
 156388 phiGeoBd[ig].resize(m_numPointBdEle[ibd]);
95
2/4
✓ Branch 1 taken 156388 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 156388 times.
✗ Branch 6 not taken.
 156388 const QuadraturePoint& quadPtRefBd = bdEle(ibd).quadratureRule().quadraturePoint(ig);
96 //std::cout << "quadPtRefBd : " ; quadPtRefBd.print(10); std::cout << std::endl;
97 /// get the shape function of the boundary element at the Boundary reference quadrature points
98
2/2
✓ Branch 1 taken 412064 times.
✓ Branch 2 taken 156388 times.
 568452 for(int iptBd=0; iptBd < m_numPointBdEle[ibd] ; iptBd++ ) {
99
3/6
✓ Branch 1 taken 412064 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 412064 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 412064 times.
✗ Branch 11 not taken.
 412064 phiGeoBd[ig](iptBd) = bdEle(ibd).geoEle().basisFunction().phi(iptBd, quadPtRefBd);
100 }
101 //std::cout << "phiGeoBd[ig=" << ig << "] = " << *phiGeoBd[ig] << std::endl;
102
103 /// compute the quad points onto the boundaries of the ref element
104 /// ("coorMap a la main" : bd points not stored)
105 /// phiGeoBd[ig] contains "weights" you use them to sum the point of the reference element
106 /// of the boundary (that have m_refCoor -1 coordinate) you obtain back quadPtRefBd.
107 /// if you use them to weigth the corresponding point of the boundary element in the
108 /// the reference space of volume element (that have m_refCoor coordinate) you obtain the same
109 /// point, but mapped from the bd reference space to the vol reference space.
110 156388 Point& theQuadPointBd = m_refQuadPointBdEle[ibd][ig];
111 156388 theQuadPointBd.x() = theQuadPointBd.y() = theQuadPointBd.z() = 0.;
112
2/2
✓ Branch 0 taken 371744 times.
✓ Branch 1 taken 156388 times.
 528132 for(int icoor=0; icoor<m_numCoor; icoor++) {
113
2/2
✓ Branch 1 taken 1041352 times.
✓ Branch 2 taken 371744 times.
 1413096 for(int iptBd=0; iptBd < m_numPointBdEle[ibd] ; iptBd++ ) {
114
2/4
✓ Branch 1 taken 1041352 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1041352 times.
✗ Branch 6 not taken.
 1041352 theQuadPointBd[icoor] += pointCoorOnBd( iptBd, icoor ) * phiGeoBd[ig](iptBd);
115 }
116 }
117 //std::cout << "currentQuadPointBd[ig=" << ig << "] = "; m_refQuadPointBdEle[ibd][ig].print(10);
118 }
119 56258 }
120
121 //--------------------
122 /// 3/ compute the indirections for accessing the quad points on the boundary
123 //--------------------
124
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_indexQuadPoint.resize(m_numBdEle+2); // indexQuadPoint[0]=0 and indexQuadPoint[1]=m_numQuadraturePoint;
125
2/2
✓ Branch 0 taken 56258 times.
✓ Branch 1 taken 16912 times.
 73170 for ( int ibd = 0; ibd <m_numBdEle; ibd++ ) {
126 56258 m_indexQuadPoint[ibd+2] = m_indexQuadPoint[ibd+1] + m_numQuadPointBdEle[ibd];
127 }
128
129 //for ( unsigned int ii = 0;ii < m_indexQuadPoint.size();ii++ ){
130 //std::cout << "m_indexQuadPoint[ii=" << ii << "] = " << m_indexQuadPoint[ii] << std::endl;}
131
132 16912 m_numQuadraturePointInternAndBd = m_indexQuadPoint.back();
133 //std::cout << "numQuadraturePointInternBd = " << m_numQuadraturePointInternAndBd << std::endl;
134
135
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_dPhiRef.resize(m_numQuadraturePointInternAndBd );
136
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_jacobian.resize(m_numQuadraturePointInternAndBd);
137
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_phiGeo.resize(m_numQuadraturePointInternAndBd);
138
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 m_dPhiGeo.resize(m_numQuadraturePointInternAndBd);
139
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 phi.resize(m_numQuadraturePointInternAndBd);
140
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 dPhi.resize(m_numQuadraturePointInternAndBd);
141
1/2
✓ Branch 1 taken 16912 times.
✗ Branch 2 not taken.
 16912 currentQuadPoint.resize(m_numQuadraturePointInternAndBd);
142
143 //--------------------
144 /// 4/ compute phi and its derivatives at the boundary quadrature points on the reference element
145 //--------------------
146
2/2
✓ Branch 0 taken 56258 times.
✓ Branch 1 taken 16912 times.
 73170 for ( int ibd = 0; ibd <m_numBdEle; ibd++ ) {
147 //std::cout << " ibd=" << ibd << std::endl;
148
2/2
✓ Branch 1 taken 156388 times.
✓ Branch 2 taken 56258 times.
 212646 for(int ilocg=0; ilocg< m_numQuadPointBdEle[ibd]; ilocg++) {
149 156388 int ig = ilocg+m_indexQuadPoint[ibd+1];
150 //std::cout << " ilocg=" << ilocg << " ig=" << ig << std::endl;
151
1/2
✓ Branch 2 taken 156388 times.
✗ Branch 3 not taken.
 156388 m_dPhiRef[ig].resize(m_numRefCoor,m_numDof);
152
1/2
✓ Branch 2 taken 156388 times.
✗ Branch 3 not taken.
 156388 m_jacobian[ig].resize(m_numCoor,m_numCoor);
153
1/2
✓ Branch 2 taken 156388 times.
✗ Branch 3 not taken.
 156388 m_phiGeo[ig].resize(m_numPoint);
154
1/2
✓ Branch 2 taken 156388 times.
✗ Branch 3 not taken.
 156388 m_dPhiGeo[ig].resize(m_numRefCoor,m_numPoint);
155
1/2
✓ Branch 2 taken 156388 times.
✗ Branch 3 not taken.
 156388 phi[ig].resize(m_numDof);
156
1/2
✓ Branch 2 taken 156388 times.
✗ Branch 3 not taken.
 156388 dPhi[ig].resize(m_numRefCoor,m_numDof);
157 //std::cout << " resize done" << std::endl << std::flush;
158
159 156388 Point& theQuadPointBd = m_refQuadPointBdEle[ibd][ilocg];
160
161 //std::cout << " setting phi/dphiref" << std::endl << std::flush;
162
2/2
✓ Branch 0 taken 667972 times.
✓ Branch 1 taken 156388 times.
 824360 for(int i=0; i<m_numDof; i++ ) {
163
2/4
✓ Branch 2 taken 667972 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 667972 times.
✗ Branch 7 not taken.
 667972 phi[ig](i) = m_refEle.basisFunction().phi(i, theQuadPointBd);
164
2/2
✓ Branch 0 taken 1710056 times.
✓ Branch 1 taken 667972 times.
 2378028 for(int icoor=0; icoor<m_numRefCoor; icoor++ ) {
165
2/4
✓ Branch 2 taken 1710056 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1710056 times.
✗ Branch 7 not taken.
 1710056 m_dPhiRef[ig](icoor,i) = m_refEle.basisFunction().dPhi(i,icoor, theQuadPointBd);
166 }
167 }
168 //std::cout << " setting phiGeo/dphiGeo" << std::endl << std::flush;
169
2/2
✓ Branch 0 taken 637732 times.
✓ Branch 1 taken 156388 times.
 794120 for(int k=0; k<m_numPoint; k++) {
170
2/4
✓ Branch 2 taken 637732 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 637732 times.
✗ Branch 7 not taken.
 637732 m_phiGeo[ig](k)= m_geoEle.basisFunction().phi(k, theQuadPointBd);
171
2/2
✓ Branch 0 taken 1619336 times.
✓ Branch 1 taken 637732 times.
 2257068 for(int icoor=0; icoor<m_numRefCoor; icoor++) {
172
2/4
✓ Branch 2 taken 1619336 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1619336 times.
✗ Branch 7 not taken.
 1619336 m_dPhiGeo[ig](icoor,k) = m_geoEle.basisFunction().dPhi(k,icoor, theQuadPointBd);
173 }
174 }
175 }
176 }
177 //std::cout << "Constructor CurrentFiniteElementWithBd done" << std::endl;
178 16912 }
179
180 /***********************************************************************************/
181 /***********************************************************************************/
182
183 6944 CurrentFiniteElementWithBd::CurrentFiniteElementWithBd(
184 const RefElement& refEle,
185 const GeoElement& geoEle,
186 const DegreeOfExactness& degOfExactness,
187 const DegreeOfExactness& degOfExactnessBd
188
2/4
✓ Branch 2 taken 6944 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6944 times.
✗ Branch 6 not taken.
 6944 ): CurrentFiniteElementWithBd(refEle,geoEle,std::vector<DegreeOfExactness>(1, degOfExactness), degOfExactnessBd)
189 {
190
191 6944 }
192
193 /***********************************************************************************/
194 /***********************************************************************************/
195
196 53968 CurrentFiniteElementWithBd::~CurrentFiniteElementWithBd()
197 {
198
2/2
✓ Branch 0 taken 56258 times.
✓ Branch 1 taken 16912 times.
 146340 for(int ibd=0; ibd<m_numBdEle; ibd++ ) {
199
1/2
✓ Branch 0 taken 56258 times.
✗ Branch 1 not taken.
 112516 delete m_bdEle[ibd];
200 }
201
1/2
✓ Branch 0 taken 16912 times.
✗ Branch 1 not taken.
 33824 delete [] m_bdEle;
202 53968 }
203
204 /***********************************************************************************/
205 /***********************************************************************************/
206
207 19667552 const CurvilinearFiniteElement & CurrentFiniteElementWithBd::bdEle(const std::size_t ibd ) const
208 {
209
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 19667552 times.
 19667552 FEL_ASSERT(ibd < static_cast<std::size_t>(m_numBdEle) );
210 19667552 const CurvilinearFiniteElement* retPtr = m_bdEle[ ibd ];
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 19667552 times.
 19667552 FEL_ASSERT(retPtr);
212 19667552 return *retPtr;
213 }
214
215 /***********************************************************************************/
216 /***********************************************************************************/
217
218 43697491 CurvilinearFiniteElement & CurrentFiniteElementWithBd::bdEle(const std::size_t ibd )
219 {
220
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43697491 times.
 43697491 FEL_ASSERT(ibd < static_cast<std::size_t>(m_numBdEle) );
221 43697491 CurvilinearFiniteElement* retPtr = m_bdEle[ ibd ];
222
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43697491 times.
 43697491 FEL_ASSERT(retPtr);
223 43697491 return *retPtr;
224 }
225
226 /***********************************************************************************/
227 /***********************************************************************************/
228
229 405472 CurvilinearFiniteElement* CurrentFiniteElementWithBd::ptrBdEle(const std::size_t ibd ) const
230 {
231
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 405472 times.
 405472 FEL_ASSERT(ibd < static_cast<std::size_t>(m_numBdEle) );
232 405472 return m_bdEle[ ibd ];
233 }
234
235 /***********************************************************************************/
236 /***********************************************************************************/
237
238 ✗ void CurrentFiniteElementWithBd::updateBdMeas()
239 {
240 ✗ FEL_ASSERT( hasPoint() ); /// you need a minimum info on the element!
241
242 ✗ for ( int ibd = 0; ibd <m_numBdEle; ibd++ ) {
243 //note: what to do of the bd elem id?? (default: local numbering in bd elem) vm 08/2011
244 ✗ CurvilinearFiniteElement* fe = m_bdEle[ibd];
245 ✗ FEL_ASSERT(fe);
246 ✗ fe->updateMeas(ibd, m_point, m_pointBdEle[ibd]);
247 }
248 }
249
250 /***********************************************************************************/
251 /***********************************************************************************/
252
253 1295565 void CurrentFiniteElementWithBd::updateBdMeasNormal()
254 {
255
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1295565 times.
 1295565 FEL_ASSERT( hasPoint() ); /// you need a minimum info on the element!
256
257
2/2
✓ Branch 0 taken 4080132 times.
✓ Branch 1 taken 1295565 times.
 5375697 for ( int ibd = 0; ibd <m_numBdEle; ibd++ ) {
258 //note: what to do of the bd elem id?? (default: local numbering in bd elem) vm 08/2011
259 4080132 CurvilinearFiniteElement* fe = m_bdEle[ibd];
260
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4080132 times.
 4080132 FEL_ASSERT(fe);
261 4080132 fe->updateMeasNormal(ibd, m_point, m_pointBdEle[ibd]);
262 }
263 1295565 }
264
265 /***********************************************************************************/
266 /***********************************************************************************/
267
268 ✗ void CurrentFiniteElementWithBd::computeCurrentQuadraturePointInternAndBd()
269 {
270 /// internal quad points:
271 ✗ CurrentFiniteElement::computeCurrentQuadraturePoint();
272 /// quad points on boundary elements: (separate treatment)
273 ✗ for (int ibd = 0; ibd <m_numBdEle; ibd++ ) {
274 ✗ for(int ilocg=0; ilocg< m_numQuadPointBdEle[ibd]; ilocg++) {
275 ✗ const std::size_t ig = ilocg+m_indexQuadPoint[ibd+1];
276 ✗ coorMap(currentQuadPoint[ig],m_refQuadPointBdEle[ibd][ilocg]);
277 }
278 }
279 }
280
281 /***********************************************************************************/
282 /***********************************************************************************/
283
284 ✗ void CurrentFiniteElementWithBd::printBd(int verbose,std::ostream& c) const
285 {
286 ✗ c << "OOoo-->CurrentFiniteElementWithBd" << std::endl;
287 ✗ c << " numBdEle: " << m_numBdEle << "\n";
288 ✗ for (int ibd = 0; ibd< m_numBdEle; ibd++ ) {
289 ✗ c << "OOoo----> Boundary Element (ibd="<< ibd << ")" << std::endl;
290 ✗ c << " Points (num=" << m_numPointBdEle[ibd] << ")";
291 ✗ if ( verbose >0) {
292 ✗ c <<" = [";
293 ✗ for(int iptBd=0; iptBd < m_numPointBdEle[ibd] ; iptBd++ ) {
294 ✗ c << m_pointBdEle[ibd][iptBd] << " ";
295 }
296 ✗ c << "]";
297 }
298 ✗ c << std::endl;
299 ✗ c << " numQuadraturePointBdEle=" << m_numQuadPointBdEle[ibd] <<std::endl;
300 ✗ if ( verbose >1) {
301 ✗ for (int ig = 0; ig < m_numQuadPointBdEle[ibd]; ig++ ) {
302 ✗ c << " ---->currentQuadPointBdEle[ig=" << ig << "] = ";
303 ✗ m_refQuadPointBdEle[ibd][ig].print(verbose,c);
304 }
305 }
306 //
307 ✗ bdEle(ibd).print(verbose,c);
308 //
309 }
310 }
311
312 /***********************************************************************************/
313 /***********************************************************************************/
314
315 // There must be an easier way to do this!! Right now I found only this one..Matteo 3/14
316 // ielbd = index of the surface element
317 // ielvol= index of the volume element
318 // dofbd = index of the local dof (->0,1,2 for a triangle)
319 // icoor = component of the variable
320 // idvar = index of the variable
321 // tested only for velocity!
322 ✗ int CurrentFiniteElementWithBd::Bd2Vol(const std::size_t iElBd, const std::size_t iElVol, const std::size_t dofBd, const std::size_t iCoor, const Dof& dof, const std::size_t idVar) const
323 {
324 felInt iDofGlob;
325 ✗ dof.loc2glob( iElBd, dofBd, idVar, iCoor, iDofGlob);
326 ✗ for ( int iVolDofL = 0; iVolDofL < m_numDof; iVolDofL++) {
327 felInt iDofGlobTry;
328 ✗ dof.loc2glob( iElVol, iVolDofL, idVar, iCoor, iDofGlobTry);
329 ✗ if( iDofGlobTry == iDofGlob)
330 ✗ return iVolDofL;
331 }
332 ✗ FEL_ERROR("Dof not found");
333 return -1;
334 }
335
336 /***********************************************************************************/
337 /***********************************************************************************/
338
339 /*
340 void CurrentFiniteElementWithBd::mapBd(int iface,double& xref3d,double& yref3d,double& zref3d,double xref2d,double yref2d){
341 GeoElement& ge=CurrentFiniteElement::m_geoEle;
342 for(int ip=0;ip<ge.numPointPerFace(); ip++){
343 for(int icoor=0;icoor<ge.numCoor();icoor++){
344 xref3d += ge.pointCoor( ge.pointOfFace(iface, ip), icoor)*;
345 }
346 }
347 }
348 */
349 }
350