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 HermiteTriangle.cpp
19 \authors D. Martin, N. Kielbasiewicz
20 \since 19 nov 2004
21 \date 6 aug 2012
22
23 \brief Implementation of xlifepp::HermiteTriangle class members and related functions
24
25 xlifepp::HermiteTriangle defines Hermite Reference Element interpolation data
26 */
27
28 #include "HermiteTriangle.hpp"
29 #include "../Interpolation.hpp"
30 #include "utils.h"
31
32 namespace xlifepp
33 {
34
35 //! HermiteTriangle constructor for Hermite reference elements
HermiteTriangle(const Interpolation * interp_p)36 HermiteTriangle::HermiteTriangle(const Interpolation* interp_p)
37 : RefTriangle(interp_p)
38 {
39 name_ += "_Her";
40 trace_p->push("HermiteTriangle::HermiteTriangle (" + name_ + ")");
41 interpolationData();
42 // local numbering of points on sides
43 sideNumbering();
44 // Reference Element on sides
45 sideRefElement();
46 maxDegree = 3;
47
48 trace_p->pop();
49 }
50
~HermiteTriangle()51 HermiteTriangle::~HermiteTriangle() {}
52
53 //! interpolationData defines Reference Element interpolation data
interpolationData()54 void HermiteTriangle::interpolationData()
55 {
56 trace_p->push("HermiteTriangle::interpolationData");
57 number_t i_n = interpolation_p->numtype;
58
59 if ( i_n == 3 ) // Hermite triangle element of degree 3
60 {
61 nbDofs_ = 10;
62 nbPts_ = 3; // 4 ?????
63 nbInternalDofs_ = 1;
64 nbDofsOnVertices_ = 9;
65 // Creating Degrees Of Freedom of reference element
66 refDofs.reserve(nbDofs_);
67 for(number_t node = 1; node <= 3; node++ )
68 {
69 refDofs.push_back(new RefDof(true, _onVertex, node, 1, 1, 0, node, 1, 0, _noProjection, "nodal value"));
70 refDofs.push_back(new RefDof(true, _onVertex, node, 1, 2, 0, node, 1, 1, _noProjection, "derivative"));
71 refDofs.push_back(new RefDof(true, _onVertex, node, 1, 3, 0, node, 1, 1, _noProjection, "derivative"));
72 }
73 refDofs.push_back(new RefDof(true, _onElement, 0, 1, 1, 0, 4, 1, 0, _noProjection, "nodal value"));
74 //hermiteRefDofs(1, 3, 2, 2); lagrangeRefDofs(10, nbDofsOnVertices_, nbDofs_ - nbInternalDofs_, 2);
75 }
76
77 // allocate storage for shape functions
78 shapeValues.resize(*this, 2);
79
80 trace_p->pop();
81 }
82
83 //! sideNumbering defines Reference Element local numbering on sides
sideNumbering()84 void HermiteTriangle::sideNumbering()
85 {
86 trace_p->push("HermiteTriangle::sideNumbering");
87 number_t i_n = interpolation_p->numtype;
88
89 number_t nbSides = geomRefElem_p->nbSides();
90 number_t more = 1;
91
92 switch (i_n)
93 {
94 case _P3: // Hermite triangle element of degree 3
95 for ( number_t side = 0; side < nbSides; side++, more += nbSides)
96 {
97 sideDofNumbers_[side].push_back(more + 3);
98 sideDofNumbers_[side].push_back(more + 4);
99 sideDofNumbers_[side].push_back(more + 5);
100 sideDofNumbers_[side].push_back(more);
101 sideDofNumbers_[side].push_back(more + 1);
102 sideDofNumbers_[side].push_back(more + 2);
103 }
104 // correction for first vertex of last edge (should be 'side+1 % nbSides')
105 sideDofNumbers_[nbSides][0] = 1;
106 sideDofNumbers_[nbSides][1] = 2;
107 sideDofNumbers_[nbSides][2] = 3;
108 default:
109 break;
110 }
111 trace_p->pop();
112 }
113
114 //! HermiteTriangle constructor for Hermite standard reference elements
115 template<number_t Pk>
HermiteStdTriangle(const Interpolation * interp_p)116 HermiteStdTriangle<Pk>::HermiteStdTriangle(const Interpolation* interp_p)
117 : HermiteTriangle(interp_p)
118 {
119 trace_p->push("HermiteTriangle::HermiteStdTriangle");
120 name_ += "_" + tostring(Pk);
121 // local coordinates of Points supporting D.o.F
122 pointCoordinates();
123
124 trace_p->pop();
125 }
126
127 template<number_t Pk>
~HermiteStdTriangle()128 HermiteStdTriangle<Pk>::~HermiteStdTriangle() {}
129
130 //! pointCoordinates defines Reference Element point coordinates
131 template<>
pointCoordinates()132 void HermiteStdTriangle<P3>::pointCoordinates()
133 {
134 trace_p->push("HermiteStdTriangle<P3>::pointCoordinates");
135 std::vector<RefDof*>::iterator it_rd = refDofs.begin();
136 // Hermite triangle element of degree 3
137 (*it_rd++)->coords(1., 0.);
138 (*it_rd++)->coords(1., 0.);
139 (*it_rd++)->coords(1., 0.);
140 (*it_rd++)->coords(0., 1.);
141 (*it_rd++)->coords(0., 1.);
142 (*it_rd++)->coords(0., 1.);
143 (*it_rd++)->coords(0., 0.);
144 (*it_rd++)->coords(0., 0.);
145 (*it_rd++)->coords(0., 0.);
146 (*it_rd)->coords(over3_, over3_);
147 trace_p->pop();
148 }
149
150 //! computeShapeValues defines Hermite Reference Element shape functions
151 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const152 void HermiteStdTriangle<P3>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
153 {
154 real_t x1 = *it_pt, x2 = *(it_pt + 1), x3 = 1. - x1 - x2;
155 real_t xx1 = x1 * (x3 - x2), xx2 = x2 * (x3 - x1), xx3 = x3 * (x2 - x1);
156 real_t Twox1 = x1 + x1, Twox2 = x2 + x2, Twox3 = x3 + x3;
157 real_t bubble = x1 * x2 * x3;
158
159 std::vector<real_t>::iterator it_w(shv.w.begin());
160 *it_w++ = x1 * x1 * (3. - Twox1) - 7 * bubble;
161 *it_w++ = -x1 * xx3;
162 *it_w++ = -x1 * xx2;
163 *it_w++ = x2 * x2 * (3. - Twox2) - 7 * bubble;
164 *it_w++ = -x2 * xx1;
165 *it_w++ = x2 * xx3;
166 *it_w++ = x3 * x3 * (3. - Twox3) - 7 * bubble;
167 *it_w++ = x3 * xx2;
168 *it_w++ = x3 * xx1;
169 *it_w++ = 27 * bubble;
170
171 if (withDeriv)
172 {
173 real_t Sevenxx1 = -7.*xx1;
174 real_t Sevenxx2 = -7.*xx2;
175 std::vector< std::vector<real_t> >::iterator it_dwdx1(shv.dw.begin()), it_dwdx2(it_dwdx1 + 1);
176 std::vector<real_t>::iterator it_dwdx1i((*it_dwdx1).begin()), it_dwdx2i((*it_dwdx2).begin());
177 *it_dwdx1i++ = 6.*x1 * (1. - x1) + Sevenxx2; *it_dwdx2i++ = Sevenxx1;
178 *it_dwdx1i++ = x1 * (Twox3 - x1) - xx2; *it_dwdx2i++ = -x1 * x1 - xx1;
179 *it_dwdx1i++ = x2 * (3.*x1 - x3); *it_dwdx2i++ = x1 * x1 - xx1;
180 *it_dwdx1i++ = Sevenxx2; *it_dwdx2i++ = 6.*x2 * (1. - x2) + Sevenxx1;
181 *it_dwdx1i++ = x2 * x2 - xx2; *it_dwdx2i++ = x1 * (3 * x2 - x3);
182 *it_dwdx1i++ = -x2 * x2 - xx2; *it_dwdx2i++ = x2 * (Twox3 - x2) - xx1;
183 *it_dwdx1i++ = -6.*x3 * (1. - x3) + Sevenxx2; *it_dwdx2i++ = -6.*x3 * (1. - x3) + Sevenxx1;
184 *it_dwdx1i++ = x2 * (x1 - 3.*x3); *it_dwdx2i++ = x3 * (x3 - Twox2) - xx1;
185 *it_dwdx1i++ = x3 * (x3 - Twox1) - xx2; *it_dwdx2i++ = x1 * (x2 - 3.*x3);
186 *it_dwdx1i = 27.*xx2; *it_dwdx2i++ = 27.*xx1;
187 }
188 /* ****
189 {
190 Twox3 = Twox3-1.;
191 shapeValues.ddw[0][0]=14.*x2-12.*x1+6.;
192 shapeValues.ddw[1][0]=-7.*Twox3;
193 shapeValues.ddw[2][0]=14.*x1;
194 shapeValues.ddw[0][1]=2.-6.*x1;
195 shapeValues.ddw[1][1]=-Twox1-Twox3;
196 shapeValues.ddw[2][1]=Twox2;
197 shapeValues.ddw[0][2]=4.*x2;
198 shapeValues.ddw[1][2]=Twox1-Twox3;
199 shapeValues.ddw[2][2]=Twox1;
200 shapeValues.ddw[0][3]=14*x2;
201 shapeValues.ddw[1][3]=-7.*Twox3;;
202 shapeValues.ddw[2][3]=14.*x1-12.*x2+6.;
203 shapeValues.ddw[0][4]=Twox2;
204 shapeValues.ddw[1][4]=2*x2-Twox3;
205 shapeValues.ddw[2][4]=4.*x1;
206 shapeValues.ddw[0][5]=Twox2;
207 shapeValues.ddw[1][5]=-Twox2-Twox3;
208 shapeValues.ddw[2][5]=2.-6.*x2;
209 shapeValues.ddw[0][6]=-6.*Twox3+14.*x2;
210 shapeValues.ddw[1][6]=-13.*Twox3;
211 shapeValues.ddw[2][6]=-6.*Twox3+14.*x1;
212 shapeValues.ddw[0][7]=4.*x2;
213 shapeValues.ddw[1][7]=x1+3.*(x2-x3);
214 shapeValues.ddw[2][7]=2.-6.*x3;
215 shapeValues.ddw[0][8]=2.-6.*x3;
216 shapeValues.ddw[1][8]=x2+3.*(x1-x3);
217 shapeValues.ddw[2][8]=4.*x1;
218 shapeValues.ddw[0][9]=-54.*x2;
219 shapeValues.ddw[1][9]=27.*Twox3;
220 shapeValues.ddw[2][9]=-54.*x1;
221 **** */
222 }
223
224 template class HermiteStdTriangle<_P3>;
225
226 } // end of namespace xlifepp
227