1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file unit_GeomElement.cpp
19 \author E. Lunéville
20 \since 05 apr 2012
21 \date 05 may 2012
22
23 Low level tests of GeomElement class and related classes.
24 Almost functionalities are checked.
25 This function may either creates a reference file storing the results (check=false)
26 or compares results to those stored in the reference file (check=true)
27 It returns reporting informations in a string
28 */
29
30 //===============================================================================
31 //library dependances
32 #include "xlife++-libs.h"
33 #include "testUtils.hpp"
34
35 //===============================================================================
36 //stl dependances
37 #include <iostream>
38 #include <fstream>
39 #include <vector>
40
41 //===============================================================================
42
43 using namespace xlifepp;
44
45 namespace unit_GeomElement {
46
unit_GeomElement(bool check)47 String unit_GeomElement(bool check)
48 {
49 String rootname = "unit_GeomElement";
50 trace_p->push(rootname);
51 std::stringstream out; //string stream receiving results
52 out.precision(testPrec);
53 verboseLevel(3);
54
55 //get interpolation structure
56 Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
57 Interpolation* interp2_p = findInterpolation(Lagrange, standard, 2, H1);
58 Number numelt = 1;
59
60 //create a segment with node 2 1
61 out << "------------------ create plain elements -------------------------------\n";
62 RefElement* ref_p = findRefElement(_segment, interp_p); //find reference element Lagrange P1
63 GeomElement S1(0, ref_p, 1, numelt);
64 S1.meshElement()->nodes[0] = new Point(0.);
65 S1.meshElement()->nodes[1] = new Point(1.);
66 S1.meshElement()->nodeNumbers[0] = 2; S1.meshElement()->nodeNumbers[1] = 1;
67 out << "Element S1 : " << S1 << std::endl;
68 numelt++;
69
70 //create a triangle with node 3 2 1
71 ref_p = findRefElement(_triangle, interp_p); //find reference element Lagrange P1
72 GeomElement T1(0, ref_p, 2, numelt);
73 T1.meshElement()->nodes[0] = new Point(0., 0.);
74 T1.meshElement()->nodes[1] = new Point(2., 0.);
75 T1.meshElement()->nodes[2] = new Point(0., 2.);
76 T1.meshElement()->nodeNumbers[0] = 3; T1.meshElement()->nodeNumbers[1] = 2; T1.meshElement()->nodeNumbers[2] = 1;
77 T1.meshElement()->vertexNumbers = T1.meshElement()->nodeNumbers;
78 out << "Element T1 : " << T1 << std::endl;
79 numelt++;
80
81 //create a P2 triangle with node 1 2 3 4 5 6
82 ref_p = findRefElement(_triangle, interp2_p); //find reference element Lagrange P1
83 GeomElement T2(0, ref_p, 2, numelt);
84 T2.meshElement()->nodes[0] = new Point(0., 0.);
85 T2.meshElement()->nodes[1] = new Point(2., 0.);
86 T2.meshElement()->nodes[2] = new Point(0., 2.);
87 T2.meshElement()->nodes[3] = new Point(1., 0.);
88 T2.meshElement()->nodes[4] = new Point(1., 1.);
89 T2.meshElement()->nodes[5] = new Point(0., 1.);
90 T2.meshElement()->nodeNumbers[0] = 1; T2.meshElement()->nodeNumbers[1] = 2; T2.meshElement()->nodeNumbers[2] = 3;
91 T2.meshElement()->nodeNumbers[3] = 4; T2.meshElement()->nodeNumbers[4] = 5; T2.meshElement()->nodeNumbers[5] = 6;
92 for(Number i = 0; i < T2.numberOfVertices(); i++)
93 { T2.meshElement()->vertexNumbers[i] = T2.meshElement()->nodeNumbers[i]; }
94 out << "Element T2 : " << T2 << std::endl;
95 numelt++;
96
97 //create side elements
98 out << "------------------- create side elements -------------------------------\n";
99 GeomElement S1_1(&S1, 1, numelt);
100 out << "Side element S1_1 : " << S1_1 << std::endl;
101 numelt++;
102 GeomElement T1_2(&T1, 2, numelt);
103 out << "Side element T1_2 : " << T1_2 << std::endl;
104 numelt++;
105 GeomElement T2_3(&T2, 3, numelt);
106 out << "Side element T2_3 : " << T2_3 << std::endl;
107 numelt++;
108
109 //test GeomElement accessors for a plain element
110 out << "------------ test accessors of GeomElement for plain elements -----------\n";
111 out << "Element T1->meshElement() : " << (*T1.meshElement()) << "\n";
112 out << "Element T1->elementDim() : " << T1.elementDim() << "\n";
113 out << "Element T1->refElement() : " << (*T1.refElement()) << "\n";
114 out << "Element T1->geomRefElement() : " << (*T1.geomRefElement()) << "\n";
115 out << "Element T1->shapeType() : " << words("shape", T1.shapeType()) << "\n";
116 out << "Element T1->numberOfVertices(), numberOfSides(), numberOfSideofSides() : "
117 << T1.numberOfVertices() << ", " << T1.numberOfSides() << ", " << T1.numberOfSideOfSides();
118 out << "\nList of vertex numbers : ";
119 for(Number i = 1; i <= T1.numberOfVertices(); i++) { out << T1.vertexNumber(i) << " "; }
120 out << "\nList of side shapes : ";
121 for(Number i = 1; i <= T1.numberOfSides(); i++) { out << i << "->" << words("shape", T1.shapeType(i)) << " "; }
122 //test specific MeshElement accessors for a plain element
123 out << "\nList of node numbers : ";
124 for(Number i = 0; i < T1.meshElement()->numberOfNodes(); i++) { out << T1.meshElement()->nodeNumbers[i] << " "; }
125 out << "\n";
126 //test GeomElement accessors for a side element
127 out << "------------ test accessors of GeomElement for side elements -----------\n";
128 out << "Element T1_2->parent() : " << (*T1_2.parent()) << "\n";
129 GeoNumPair gnp = T1_2.parentSide(0);
130 out << "Element T1_2->parentSide(0).second : " << gnp.second;
131 out << " Element T1_2->elementDim() : " << T1_2.elementDim() << "\n";
132 out << "Element T1_2->refElement() : " << (*T1_2.refElement()) << "\n";
133 out << "Element T1_2->geomRefElement() : " << (*T1_2.geomRefElement()) << "\n";
134 out << "Element T1_2->shapeType() : " << words("shape", T1_2.shapeType()) << "\n";
135 out << "Element T1_2->numberOfVertices(), numberOfSides(), numberOfSideofSides() : "
136 << T1_2.numberOfVertices() << ", " << T1_2.numberOfSides() << ", " << T1_2.numberOfSideOfSides();
137 out << "\nList of vertex numbers : ";
138 for(Number i = 1; i <= T1_2.numberOfVertices(); i++) { out << T1_2.vertexNumber(i) << " "; }
139 out << "\n";
140
141 out << "Element T2_3->parent() : " << (*T2_3.parent()) << "\n";
142 gnp = T2_3.parentSide(0);
143 out << "Element T2_3->parentSide(0).second : " << gnp.second;
144 out << " Element T2_3->elementDim() : " << T2_3.elementDim() << "\n";
145 out << "Element T2_3->refElement() : " << (*T2_3.refElement()) << "\n";
146 out << "Element T2_3->geomRefElement() : " << (*T2_3.geomRefElement()) << "\n";
147 out << "Element T2_3->shapeType() : " << words("shape", T2_3.shapeType()) << "\n";
148 out << "Element T2_3->numberOfVertices(), numberOfSides(), numberOfSideofSides() : "
149 << T2_3.numberOfVertices() << ", " << T2_3.numberOfSides() << ", " << T2_3.numberOfSideOfSides() << "\n";
150 out << "List of vertex numbers : ";
151 for(Number i = 1; i <= T2_3.numberOfVertices(); i++) { out << T2_3.vertexNumber(i) << " "; }
152 //test construction of MeshElement for a side element
153 out << "\n------------ test buildSideElement for side elements P1 and P2 -----------\n";
154 T1_2.buildSideMeshElement();
155 out << "Element T1_2 : " << T1_2 << "\n";
156 T2_3.buildSideMeshElement();
157 out << "Element T2_3 : " << T2_3 << "\n";
158
159 //test of MeshElement functions
160 MeshElement* melt = T1.meshElement();
161 out << "\n------------ test MeshElement functions P1 -----------\n";
162 out << "melt->index()=" << melt->index() << " melt->spaceDim()=" << melt->spaceDim();
163 out << " melt->shapeType()= " << words("shape", melt->shapeType()) << " melt->numberOfNodes()= " << melt->numberOfNodes();
164 out << " melt->numberOfSides()= " << melt->numberOfSides() << " melt->numberOfSideOfSides() : " << melt->numberOfSideOfSides() << "\n";
165 out << "elt.refElement()=" << (*melt->refElement()) << "\n";
166 out << "melt->geomRefElement() : " << (*melt->geomRefElement()) << "\n";
167 out << "List of vertex numbers : ";
168 for(Number i = 1; i <= melt->numberOfVertices(); i++) { out << melt->vertexNumber(i) << " "; }
169 for(Number s = 1; s <= melt->numberOfSides(); s++)
170 {
171 out << "\nside " << s << " : melt->shapeType(s) = " << words("shape", melt->shapeType(s));
172 out << " vertex numbers : ";
173 for(Number i = 0; i < melt->verticesNumbers(s).size(); i++) { out << melt->verticesNumbers(s)[i] << " "; }
174 out << "\nmelt->refElement(s) : " << (*melt->refElement(s)) << "\n";
175 out << "melt->geomRefElement()(s) : " << (*melt->geomRefElement(s));
176 }
177 melt->computeMeasures();
178 out << "\nmelt->computeMeasures() =";
179 for(Number i = 0; i <= melt->numberOfSides(); i++) { out << " " << melt->measure(i); }
180 melt->computeOrientation();
181 out << "\nmelt->computeOrientation() = " << melt->orientation;
182
183 melt = T2.meshElement();
184 out << "\n------------ test MeshElement functions P2 -----------\n";
185 out << "melt->index()=" << melt->index() << " melt->spaceDim()=" << melt->spaceDim();
186 out << " melt->shapeType()= " << words("shape", melt->shapeType()) << " melt->numberOfNodes()= " << melt->numberOfNodes();
187 out << " melt->numberOfSides()= " << melt->numberOfSides() << " melt->numberOfSideOfSides() : " << melt->numberOfSideOfSides() << "\n";
188 out << "elt.refElement()=" << (*melt->refElement()) << "\n";
189 out << "melt->geomRefElement() : " << (*melt->geomRefElement()) << "\n";
190 out << "List of vertex numbers : ";
191 for(Number i = 1; i <= melt->numberOfVertices(); i++) { out << melt->vertexNumber(i) << " "; }
192 for(Number s = 1; s <= melt->numberOfSides(); s++)
193 {
194 out << "\nside " << s << " : melt->shapeType(s) = " << words("shape", melt->shapeType(s));
195 out << " vertex numbers : ";
196 for(Number i = 0; i < melt->verticesNumbers(s).size(); i++) { out << melt->verticesNumbers(s)[i] << " "; }
197 out << "\nmelt->refElement(s) : " << (*melt->refElement(s)) << "\n";
198 out << "melt->geomRefElement()(s) : " << (*melt->geomRefElement(s));
199 }
200 melt->computeMeasures();
201 out << "\nmelt->computeMeasures() =";
202 for(Number i = 0; i <= melt->numberOfSides(); i++) { out << " " << melt->measure(i); }
203 melt->computeOrientation();
204 out << "\nmelt->computeOrientation() = " << melt->orientation;
205
206 out << "\n------------ test GeomMapData functions P1-----------\n";
207 GeomMapData mapdata(T1.meshElement());
208 out << "mapdata.geomMap(Point(0.,0.)) = " << mapdata.geomMap(Point(0., 0.));
209 out << "\nmapdata.geomMap(Point(1.,0.)) = " << mapdata.geomMap(Point(1., 0.));
210 out << "\nmapdata.geomMap(Point(0.,1.)) = " << mapdata.geomMap(Point(0., 1.));
211 out << "\nmapdata.geomMap(Point(1./3,1./3)) = " << mapdata.geomMap(Point(1. / 3, 1. / 3));
212 mapdata.computeJacobianMatrix(Point(0., 0.));
213 out << "\nmapdata.computeJacobianMatrix(Point(0,0)) = " << mapdata.jacobianMatrix;
214 out << " determinant = " << mapdata.computeJacobianDeterminant();
215 mapdata.invertJacobianMatrix();
216 out << "\nmapdata.invertJacobianMatrix() = " << mapdata.inverseJacobianMatrix;
217 mapdata.computeJacobianMatrix(Point(0., 0.), 1);
218 out << "\nmapdata.computeJacobianMatrix(Point(0,0),1) = " << mapdata.jacobianMatrix;
219
220 out << "\n------------ test GeomMapData functions P2-----------\n";
221 mapdata = GeomMapData(T2.meshElement());
222 out << "mapdata.geomMap(Point(0.,0.)) = " << mapdata.geomMap(Point(0., 0.));
223 out << "\nmapdata.geomMap(Point(1.,0.)) = " << mapdata.geomMap(Point(1., 0.));
224 out << "\nmapdata.geomMap(Point(0.,1.)) = " << mapdata.geomMap(Point(0., 1.));
225 out << "\nmapdata.geomMap(Point(1./3,1./3)) = " << mapdata.geomMap(Point(1. / 3, 1. / 3));
226 mapdata.computeJacobianMatrix(Point(0., 0.));
227 out << "\nmapdata.computeJacobianMatrix(Point(0,0)) = " << mapdata.jacobianMatrix;
228 out << " determinant = " << mapdata.computeJacobianDeterminant();
229 mapdata.invertJacobianMatrix();
230 out << "\nmapdata.invertJacobianMatrix() = " << mapdata.inverseJacobianMatrix;
231 mapdata.computeJacobianMatrix(Point(0., 0.), 1);
232 out << "\nmapdata.computeJacobianMatrix(Point(0,0),1) = " << mapdata.jacobianMatrix;
233
234 //------------------------------------------------------------------------------------
235 // save results in a file or compare results with some references value in a file
236 //------------------------------------------------------------------------------------
237 trace_p->pop();
238 if (check) { return diffResFile(out, rootname); }
239 else { return saveResToFile(out, rootname); }
240
241 }
242
243 }
244