Intrepid
Intrepid_HCURL_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 
50 namespace 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  const int littleN = n*(n+1); // dim of (P_{n-1})^2 -- smaller space
67  const int bigN = (n+1)*(n+2); // dim of (P_{n})^2 -- larger space
68  const int scalarSmallestN = (n-1)*n / 2;
69  const int scalarLittleN = littleN/2;
70  const int scalarBigN = bigN/2;
71 
72  // first, need to project the basis for Nedelec space onto the
73  // orthogonal basis of degree n
74  // get coefficients of PkHx
75 
76  Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, N);
77 
78  // basis for the space is
79  // { (phi_i,0) }_{i=0}^{scalarLittleN-1} ,
80  // { (0,phi_i) }_{i=0}^{scalarLittleN-1} ,
81  // { (x,y) \times phi_i}_{i=scalarLittleN}^{scalarBigN-1}
82  // { (x,y) \times phi = (y phi , -x \phi)
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) \times 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,1)
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,0)
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  if (pointType == POINTTYPE_WARPBLEND) {
135  CubatureDirectLineGauss<Scalar> edgeRule( 2*n - 1 );
136  FieldContainer<Scalar> edgeCubWts( n );
137  edgeRule.getCubature( linePts , edgeCubWts );
138  }
139  else if (pointType == POINTTYPE_EQUISPACED ) {
140  shards::CellTopology linetop(shards::getCellTopologyData<shards::Line<2> >() );
141 
142  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( linePts ,
143  linetop ,
144  n+1 , 1 ,
145  POINTTYPE_EQUISPACED );
146  }
147 
148 
149  FieldContainer<Scalar> edgePts( n , 2 );
150  FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , n );
151  FieldContainer<Scalar> edgeTan(2);
152 
153  for (int i=0;i<3;i++) { // loop over edges
155  i ,
156  this->basisCellTopology_ );
157  /* multiply by 2.0 to account for a Jacobian in Pavel's definition */
158  for (int j=0;j<2;j++) {
159  edgeTan(j) *= 2.0;
160  }
161 
163  linePts ,
164  1 ,
165  i ,
166  this->basisCellTopology_ );
167 
168  Phis_.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE );
169 
170  // loop over points (rows of V2)
171  for (int j=0;j<n;j++) {
172  // loop over orthonormal basis functions (columns of V2)
173  for (int k=0;k<scalarBigN;k++) {
174  V2(n*i+j,k) = edgeTan(0) * phisAtEdgePoints(k,j);
175  V2(n*i+j,k+scalarBigN) = edgeTan(1) * phisAtEdgePoints(k,j);
176  }
177  }
178  }
179 
180  // remaining nodes are x- and y- components at internal points, if n > 1
181  // this code is exactly the same as it is for HDIV
182 
183  const int numInternalPoints = PointTools::getLatticeSize( this->getBaseCellTopology() ,
184  n + 1 ,
185  1 );
186 
187  if (numInternalPoints > 0) {
188  FieldContainer<Scalar> internalPoints( numInternalPoints , 2 );
189  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( internalPoints ,
190  this->getBaseCellTopology() ,
191  n + 1 ,
192  1 ,
193  pointType );
194 
195  FieldContainer<Scalar> phisAtInternalPoints( scalarBigN , numInternalPoints );
196  Phis_.getValues( phisAtInternalPoints , internalPoints , OPERATOR_VALUE );
197 
198  // copy values into right positions of V2
199  for (int i=0;i<numInternalPoints;i++) {
200  for (int j=0;j<scalarBigN;j++) {
201  // x component
202  V2(3*n+i,j) = phisAtInternalPoints(j,i);
203  // y component
204  V2(3*n+numInternalPoints+i,scalarBigN+j) = phisAtInternalPoints(j,i);
205  }
206  }
207  }
208 // std::cout << "Nodes on big basis\n";
209 // std::cout << V2 << "\n";
210 // std::cout << "End nodes\n";
211 
212  Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
213 
214  // multiply V2 * V1 --> V
215  Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , V1 , 0.0 );
216 
217 // std::cout << "Vandermonde:\n";
218 // std::cout << Vsdm << "\n";
219 // std::cout << "End Vandermonde\n";
220 
221  Teuchos::SerialDenseSolver<int,Scalar> solver;
222  solver.setMatrix( rcp( &Vsdm , false ) );
223  solver.invert( );
224 
225  Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
226  Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V1 , Vsdm , 0.0 );
227 
228  // std::cout << Csdm << "\n";
229 
230  for (int i=0;i<bigN;i++) {
231  for (int j=0;j<N;j++) {
232  coeffs_(i,j) = Csdm(i,j);
233  }
234  }
235  }
236 
237  template<class Scalar, class ArrayScalar>
239 
240  // Basis-dependent initializations
241  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
242  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
243  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
244  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
245 
246  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
247 
248  int *tags = new int[ tagSize * this->getCardinality() ];
249  int *tag_cur = tags;
250  const int degree = this->getDegree();
251 
252  // there are degree internal dofs on each edge -- normals. Let's do them
253  for (int ed=0;ed<3;ed++) {
254  for (int i=0;i<degree;i++) {
255  tag_cur[0] = 1; tag_cur[1] = ed; tag_cur[2] = i; tag_cur[3] = degree;
256  tag_cur += tagSize;
257  }
258  }
259 
260  // end edge dofs
261 
262  // the rest of the dofs are internal
263  const int numFaceDof = (degree-1)*degree;
264  int faceDofCur = 0;
265  for (int i=3*degree;i<degree*(degree+2);i++) {
266  tag_cur[0] = 2; tag_cur[1] = 0; tag_cur[2] = faceDofCur; tag_cur[3] = numFaceDof;
267  tag_cur += tagSize;
268  faceDofCur++;
269  }
270 
271 
272  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
273  this -> ordinalToTag_,
274  tags,
275  this -> basisCardinality_,
276  tagSize,
277  posScDim,
278  posScOrd,
279  posDfOrd);
280 
281  delete []tags;
282 
283  }
284 
285 
286 
287  template<class Scalar, class ArrayScalar>
289  const ArrayScalar & inputPoints,
290  const EOperator operatorType) const {
291 
292  // Verify arguments
293 #ifdef HAVE_INTREPID_DEBUG
294  Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
295  inputPoints,
296  operatorType,
297  this -> getBaseCellTopology(),
298  this -> getCardinality() );
299 #endif
300  const int numPts = inputPoints.dimension(0);
301  const int deg = this -> getDegree();
302  const int scalarBigN = (deg+1)*(deg+2)/2;
303 
304  try {
305  switch (operatorType) {
306  case OPERATOR_VALUE:
307  {
308  FieldContainer<Scalar> phisCur( scalarBigN , numPts );
309  Phis_.getValues( phisCur , inputPoints , OPERATOR_VALUE );
310 
311  for (int i=0;i<outputValues.dimension(0);i++) { // RT bf
312  for (int j=0;j<outputValues.dimension(1);j++) { // point
313  outputValues(i,j,0) = 0.0;
314  outputValues(i,j,1) = 0.0;
315  for (int k=0;k<scalarBigN;k++) { // Dubiner bf
316  outputValues(i,j,0) += coeffs_(k,i) * phisCur(k,j);
317  outputValues(i,j,1) += coeffs_(k+scalarBigN,i) * phisCur(k,j);
318  }
319  }
320  }
321  }
322  break;
323  case OPERATOR_CURL:
324  {
325  FieldContainer<Scalar> phisCur( scalarBigN , numPts , 2 );
326  Phis_.getValues( phisCur , inputPoints , OPERATOR_GRAD );
327  for (int i=0;i<outputValues.dimension(0);i++) { // bf loop
328  for (int j=0;j<outputValues.dimension(1);j++) { // point loop
329  // dy of x component
330  outputValues(i,j) = 0.0;
331  for (int k=0;k<scalarBigN;k++) {
332  outputValues(i,j) -= coeffs_(k,i) * phisCur(k,j,1);
333  }
334  // -dx of y component
335  for (int k=0;k<scalarBigN;k++) {
336  outputValues(i,j) += coeffs_(k+scalarBigN,i) * phisCur(k,j,0);
337  }
338  }
339  }
340  }
341  break;
342  default:
343  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
344  ">>> ERROR (Basis_HCURL_TRI_In_FEM): Operator type not implemented");
345  break;
346  }
347  }
348  catch (std::invalid_argument &exception){
349  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
350  ">>> ERROR (Basis_HCURL_TRI_In_FEM): Operator type not implemented");
351  }
352 
353  }
354 
355 
356 
357  template<class Scalar, class ArrayScalar>
359  const ArrayScalar & inputPoints,
360  const ArrayScalar & cellVertices,
361  const EOperator operatorType) const {
362  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
363  ">>> ERROR (Basis_HCURL_TRI_In_FEM): FEM Basis calling an FVD member function");
364  }
365 
366 
367 }// namespace Intrepid
368 
369 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
370 #ifdef __GNUC__
371 #warning "The Intrepid package is deprecated"
372 #endif
373 #endif
374 
static void getReferenceEdgeTangent(ArrayEdgeTangent &refEdgeTangent, const int &edgeOrd, const shards::CellTopology &parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
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.
Basis_HGRAD_TRI_Cn_FEM_ORTH< Scalar, FieldContainer< Scalar > > Phis_
Orthogonal basis of ofder n, in terms of which the H(curl) basis functions are expressed.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Defines Gauss integration rules on a line.
EBasis basisType_
Type of the basis.
Basis_HCURL_TRI_In_FEM(const int n, const EPointType pointType)
Constructor.
bool basisTagsAreSet_
&quot;true&quot; if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual int getNumPoints() const
Returns the number of cubature points.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
FieldContainer< Scalar > coeffs_
Array holding the expansion coefficients of the nodal basis in terms of Phis_.
Defines direct integration rules on a triangle.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
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...
virtual void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.