Intrepid
Intrepid_HDIV_TRI_In_FEMDef.hpp
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
50namespace Intrepid {
51
52 template<class Scalar, class ArrayScalar>
54 const EPointType pointType ):
55 Phis( n ),
56 coeffs( (n+1)*(n+2) , n*(n+2) )
57 {
58 const int N = n*(n+2);
59 this -> basisCardinality_ = N;
60 this -> basisDegree_ = n;
61 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
62 this -> basisType_ = BASIS_FEM_FIAT;
63 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
64 this -> basisTagsAreSet_ = false;
65
66
67 const int littleN = n*(n+1); // dim of (P_{n-1})^2 -- smaller space
68 const int bigN = (n+1)*(n+2); // dim of (P_{n})^2 -- larger space
69 const int scalarSmallestN = (n-1)*n / 2;
70 const int scalarLittleN = littleN/2;
71 const int scalarBigN = bigN/2;
72
73 // first, need to project the basis for RT space onto the
74 // orthogonal basis of degree n
75 // get coefficients of PkHx
76
77 Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, N);
78
79 // basis for the space is
80 // { (phi_i,0) }_{i=0}^{scalarLittleN-1} ,
81 // { (0,phi_i) }_{i=0}^{scalarLittleN-1} ,
82 // { (x,y) . phi_i}_{i=scalarSmallestN}^{scalarLittleN-1}
83 // columns of V1 are expansion of this basis in terms of the basis
84 // for P_{n}^2
85
86 // these two loops get the first two sets of basis functions
87 for (int i=0;i<scalarLittleN;i++) {
88 V1(i,i) = 1.0;
89 V1(scalarBigN+i,scalarLittleN+i) = 1.0;
90 }
91
92 // now I need to integrate { (x,y) phi } against the big basis
93 // first, get a cubature rule.
95 FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 2 );
96 FieldContainer<Scalar> cubWeights( myCub.getNumPoints() );
97 myCub.getCubature( cubPoints , cubWeights );
98
99 // tabulate the scalar orthonormal basis at cubature points
100 FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() );
101 Phis.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE );
102
103 // now do the integration
104 for (int i=0;i<n;i++) {
105 for (int j=0;j<scalarBigN;j++) { // int (x,y) phi_i \cdot (phi_j,0)
106 V1(j,littleN+i) = 0.0;
107 for (int k=0;k<myCub.getNumPoints();k++) {
108 V1(j,littleN+i) +=
109 cubWeights(k) * cubPoints(k,0)
110 * phisAtCubPoints(scalarSmallestN+i,k)
111 * phisAtCubPoints(j,k);
112 }
113 }
114 for (int j=0;j<scalarBigN;j++) { // int (x,y) phi_i \cdot (0,phi_j)
115 V1(j+scalarBigN,littleN+i) = 0.0;
116 for (int k=0;k<myCub.getNumPoints();k++) {
117 V1(j+scalarBigN,littleN+i) +=
118 cubWeights(k) * cubPoints(k,1)
119 * phisAtCubPoints(scalarSmallestN+i,k)
120 * phisAtCubPoints(j,k);
121 }
122 }
123 }
124
125 //std::cout << V1 << "\n";
126
127
128 // next, apply the RT nodes (rows) to the basis for (P_n)^2 (columns)
129 Teuchos::SerialDenseMatrix<int,Scalar> V2(N , bigN);
130
131 // first 3 * degree nodes are normals at each edge
132 // get the points on the line
133 FieldContainer<Scalar> linePts( n , 1 );
134
135 // change by Nate Roberts 8/25/16: use getLattice() for warp blend points, too
136 // (under previous approach--which used Gauss cubature points on the line--the
137 // resulting RT basis would be numerically linearly dependent for orders >= 5.)
138 shards::CellTopology linetop(shards::getCellTopologyData<shards::Line<2> >() );
139
140 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( linePts ,
141 linetop ,
142 n+1 , 1 ,
143 pointType );
144 // holds the image of the line points
145 FieldContainer<Scalar> edgePts( n , 2 );
146 FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , n );
147
148 // these are scaled by the appropriate edge lengths.
149 const Scalar nx[] = {0.0,1.0,-1.0};
150 const Scalar ny[] = {-1.0,1.0,0.0};
151
152 for (int i=0;i<3;i++) { // loop over edges
154 linePts ,
155 1 ,
156 i ,
157 this->basisCellTopology_ );
158
159 Phis.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE );
160
161 // loop over points (rows of V2)
162 for (int j=0;j<n;j++) {
163 // loop over orthonormal basis functions (columns of V2)
164 for (int k=0;k<scalarBigN;k++) {
165 V2(n*i+j,k) = nx[i] * phisAtEdgePoints(k,j);
166 V2(n*i+j,k+scalarBigN) = ny[i] * phisAtEdgePoints(k,j);
167 }
168 }
169 }
170
171 // next map the points to each edge
172
173
174 // remaining nodes are divided into two pieces: point value of x
175 // components and point values of y components. These are
176 // evaluated at the interior of a lattice of degree + 1, For then
177 // the degree == 1 space corresponds classicaly to RT0 and so gets
178 // no internal nodes, and degree == 2 corresponds to RT1 and needs
179 // one internal node per vector component.
180 const int numInternalPoints = PointTools::getLatticeSize( this->getBaseCellTopology() ,
181 n + 1 ,
182 1 );
183
184 if (numInternalPoints > 0) {
185 FieldContainer<Scalar> internalPoints( numInternalPoints , 2 );
186 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( internalPoints ,
187 this->getBaseCellTopology() ,
188 n + 1 ,
189 1 ,
190 pointType );
191
192 FieldContainer<Scalar> phisAtInternalPoints( scalarBigN , numInternalPoints );
193 Phis.getValues( phisAtInternalPoints , internalPoints , OPERATOR_VALUE );
194
195 // copy values into right positions of V2
196 for (int i=0;i<numInternalPoints;i++) {
197 for (int j=0;j<scalarBigN;j++) {
198 // x component
199 V2(3*n+i,j) = phisAtInternalPoints(j,i);
200 // y component
201 V2(3*n+numInternalPoints+i,scalarBigN+j) = phisAtInternalPoints(j,i);
202 }
203 }
204 }
205// std::cout << "Nodes on big basis\n";
206// std::cout << V2 << "\n";
207// std::cout << "End nodes\n";
208
209 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
210
211 // multiply V2 * V1 --> V
212 Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , V1 , 0.0 );
213
214// std::cout << "Vandermonde:\n";
215// std::cout << Vsdm << "\n";
216// std::cout << "End Vandermonde\n";
217
218 Teuchos::SerialDenseSolver<int,Scalar> solver;
219 solver.setMatrix( rcp( &Vsdm , false ) );
220 solver.invert( );
221
222 Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
223 Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V1 , Vsdm , 0.0 );
224
225 // std::cout << Csdm << "\n";
226
227 for (int i=0;i<bigN;i++) {
228 for (int j=0;j<N;j++) {
229 coeffs(i,j) = Csdm(i,j);
230 }
231 }
232 }
233
234 template<class Scalar, class ArrayScalar>
236
237 // Basis-dependent initializations
238 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
239 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
240 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
241 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
242
243 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
244
245 int *tags = new int[ tagSize * this->getCardinality() ];
246 int *tag_cur = tags;
247 const int degree = this->getDegree();
248
249 // there are degree internal dofs on each edge -- normals. Let's do them
250 for (int ed=0;ed<3;ed++) {
251 for (int i=0;i<degree;i++) {
252 tag_cur[0] = 1; tag_cur[1] = ed; tag_cur[2] = i; tag_cur[3] = degree;
253 tag_cur += tagSize;
254 }
255 }
256
257 // end edge dofs
258
259 // the rest of the dofs are internal
260 const int numFaceDof = (degree-1)*degree;
261 int faceDofCur = 0;
262 for (int i=3*degree;i<3*degree+numFaceDof;i++) {
263 tag_cur[0] = 2; tag_cur[1] = 0; tag_cur[2] = faceDofCur; tag_cur[3] = numFaceDof;
264 tag_cur += tagSize;
265 faceDofCur++;
266 }
267
268
269 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
270 this -> ordinalToTag_,
271 tags,
272 this -> basisCardinality_,
273 tagSize,
274 posScDim,
275 posScOrd,
276 posDfOrd);
277
278 delete []tags;
279
280 }
281
282
283
284 template<class Scalar, class ArrayScalar>
286 const ArrayScalar & inputPoints,
287 const EOperator operatorType) const {
288
289 // Verify arguments
290#ifdef HAVE_INTREPID_DEBUG
291 Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
292 inputPoints,
293 operatorType,
294 this -> getBaseCellTopology(),
295 this -> getCardinality() );
296#endif
297 const int numPts = inputPoints.dimension(0);
298 const int deg = this -> getDegree();
299 const int scalarBigN = (deg+1)*(deg+2)/2;
300
301 try {
302 switch (operatorType) {
303 case OPERATOR_VALUE:
304 {
305 FieldContainer<Scalar> phisCur( scalarBigN , numPts );
306 Phis.getValues( phisCur , inputPoints , OPERATOR_VALUE );
307
308 for (int i=0;i<outputValues.dimension(0);i++) { // RT bf
309 for (int j=0;j<outputValues.dimension(1);j++) { // point
310 outputValues(i,j,0) = 0.0;
311 outputValues(i,j,1) = 0.0;
312 for (int k=0;k<scalarBigN;k++) { // Dubiner bf
313 outputValues(i,j,0) += coeffs(k,i) * phisCur(k,j);
314 outputValues(i,j,1) += coeffs(k+scalarBigN,i) * phisCur(k,j);
315 }
316 }
317 }
318 }
319 break;
320 case OPERATOR_DIV:
321 {
322 FieldContainer<Scalar> phisCur( scalarBigN , numPts , 2 );
323 Phis.getValues( phisCur , inputPoints , OPERATOR_GRAD );
324 for (int i=0;i<outputValues.dimension(0);i++) { // bf loop
325 for (int j=0;j<outputValues.dimension(1);j++) { // point loop
326 // dx of x component
327 outputValues(i,j) = 0.0;
328 for (int k=0;k<scalarBigN;k++) {
329 outputValues(i,j) += coeffs(k,i) * phisCur(k,j,0);
330 }
331 // dy of y component
332 for (int k=0;k<scalarBigN;k++) {
333 outputValues(i,j) += coeffs(k+scalarBigN,i) * phisCur(k,j,1);
334 }
335 }
336 }
337 }
338 break;
339 default:
340 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
341 ">>> ERROR (Basis_HDIV_TRI_In_FEM): Operator type not implemented");
342 break;
343 }
344 }
345 catch (std::invalid_argument &exception){
346 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
347 ">>> ERROR (Basis_HDIV_TRI_In_FEM): Operator type not implemented");
348 }
349
350 }
351
352
353
354 template<class Scalar, class ArrayScalar>
356 const ArrayScalar & inputPoints,
357 const ArrayScalar & cellVertices,
358 const EOperator operatorType) const {
359 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
360 ">>> ERROR (Basis_HDIV_TRI_In_FEM): FEM Basis calling an FVD member function");
361 }
362
363
364}// namespace Intrepid
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
Basis_HDIV_TRI_In_FEM(const int n, const EPointType pointType)
Constructor.
FieldContainer< Scalar > coeffs
expansion coefficients of the nodal basis in terms of the orthgonal one
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
Basis_HGRAD_TRI_Cn_FEM_ORTH< Scalar, FieldContainer< Scalar > > Phis
Orthogonal basis out of which the nodal basis is constructed.
virtual void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
EBasis basisType_
Type of the basis.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos....
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
Defines direct integration rules on a triangle.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
virtual int getNumPoints() const
Returns the number of cubature points.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
static int getLatticeSize(const shards::CellTopology &cellType, const int order, const int offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...