Intrepid
Intrepid_HGRAD_TRI_Cn_FEMDef.hpp
1 #ifndef INTREPID_HGRAD_TRI_CN_FEMDEF_HPP
2 #define INTREPID_HGRAD_TRI_CN_FEMDEF_HPP
3 // @HEADER
4 // ************************************************************************
5 //
6 // Intrepid Package
7 // Copyright (2007) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
40 // Denis Ridzal (dridzal@sandia.gov), or
41 // Kara Peterson (kjpeter@sandia.gov)
42 //
43 // ************************************************************************
44 // @HEADER
45 
51 namespace Intrepid {
52 
53  template<class Scalar, class ArrayScalar>
55  const EPointType pointType ):
56  Phis( n ),
57  V((n+1)*(n+2)/2,(n+1)*(n+2)/2),
58  Vinv((n+1)*(n+2)/2,(n+1)*(n+2)/2),
59  latticePts( (n+1)*(n+2)/2 , 2 )
60  {
61  TEUCHOS_TEST_FOR_EXCEPTION( n <= 0, std::invalid_argument, "polynomial order must be >= 1");
62 
63  const int N = (n+1)*(n+2)/2;
64  this -> basisCardinality_ = N;
65  this -> basisDegree_ = n;
66  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
67  this -> basisType_ = BASIS_FEM_FIAT;
68  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
69  this -> basisTagsAreSet_ = false;
70 
71  // construct lattice
72 
73  shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() );
74 
75  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts ,
76  myTri_3 ,
77  n ,
78  0 ,
79  pointType );
80 
81 
82  // form Vandermonde matrix. Actually, this is the transpose of the VDM,
83  // so we transpose on copy below.
84 
85  Phis.getValues( V , latticePts , OPERATOR_VALUE );
86 
87  // now I need to copy V into a Teuchos array to do the inversion
88  Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
89  for (int i=0;i<N;i++) {
90  for (int j=0;j<N;j++) {
91  Vsdm(i,j) = V(i,j);
92  }
93  }
94 
95  // invert the matrix
96  Teuchos::SerialDenseSolver<int,Scalar> solver;
97  solver.setMatrix( rcp( &Vsdm , false ) );
98  solver.invert( );
99 
100  // now I need to copy the inverse into Vinv
101  for (int i=0;i<N;i++) {
102  for (int j=0;j<N;j++) {
103  Vinv(i,j) = Vsdm(j,i);
104  }
105  }
106 
107  }
108 
109 
110  template<class Scalar, class ArrayScalar>
111  void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::initializeTags() {
112 
113  // Basis-dependent initializations
114  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
115  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
116  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
117  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
118 
119  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
120 
121  int *tags = new int[ tagSize * this->getCardinality() ];
122  int *tag_cur = tags;
123  const int degree = this->getDegree();
124 
125  // BEGIN DOF ALONG BOTTOM EDGE
126 
127  // the first dof is on vertex 0
128  tag_cur[0] = 0; tag_cur[1] = 0; tag_cur[2] = 0; tag_cur[3] = 1;
129  tag_cur += tagSize;
130 
131  // next degree-1 dof are on edge 0
132  for (int i=1;i<degree;i++) {
133  tag_cur[0] = 1; tag_cur[1] = 0; tag_cur[2] = i-1; tag_cur[3] = degree-1;
134  tag_cur += tagSize;
135  }
136 
137  // last dof is on vertex 1
138  tag_cur[0] = 0; tag_cur[1] = 1; tag_cur[2] = 0; tag_cur[3] = 1;
139  tag_cur += tagSize;
140 
141  // END DOF ALONG BOTTOM EDGE
142 
143  int num_internal_dof = PointTools::getLatticeSize( this->getBaseCellTopology() ,
144  this->getDegree() ,
145  1 );
146 
147  int internal_dof_cur = 0;
148 
149  // BEGIN DOF ALONG INTERNAL HORIZONTAL LINES
150  for (int i=1;i<degree;i++) {
151  // first dof along line is on edge #2
152  tag_cur[0] = 1; tag_cur[1] = 2; tag_cur[2] = i-1; tag_cur[3] = degree-1;
153  tag_cur += tagSize;
154 
155  // next dof are internal
156  for (int j=1;j<degree-i;j++) {
157  tag_cur[0] = 2; tag_cur[1] = 0; tag_cur[2] = internal_dof_cur; tag_cur[3] = num_internal_dof;
158  internal_dof_cur++;
159  tag_cur += tagSize;
160  }
161 
162  // last dof along line is on edge 1
163  tag_cur[0] = 1; tag_cur[1] = 1; tag_cur[2] = i-1; tag_cur[3] = degree-1;
164  tag_cur += tagSize;
165 
166  }
167  // END DOF ALONG INTERNAL HORIZONTAL LINES
168 
169  // LAST DOF IS ON VERTEX 2
170  tag_cur[0] = 0; tag_cur[1] = 2; tag_cur[2] = 0; tag_cur[3] = 1;
171  // END LAST DOF
172 
173 
174  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
175  this -> ordinalToTag_,
176  tags,
177  this -> basisCardinality_,
178  tagSize,
179  posScDim,
180  posScOrd,
181  posDfOrd);
182 
183  delete []tags;
184 
185  }
186 
187 
188 
189  template<class Scalar, class ArrayScalar>
190  void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues,
191  const ArrayScalar & inputPoints,
192  const EOperator operatorType) const {
193 
194  // Verify arguments
195 #ifdef HAVE_INTREPID_DEBUG
196  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
197  inputPoints,
198  operatorType,
199  this -> getBaseCellTopology(),
200  this -> getCardinality() );
201 #endif
202  const int numPts = inputPoints.dimension(0);
203  const int numBf = this->getCardinality();
204 
205  try {
206  switch (operatorType) {
207  case OPERATOR_VALUE:
208  {
209  FieldContainer<Scalar> phisCur( numBf , numPts );
210  Phis.getValues( phisCur , inputPoints , operatorType );
211  for (int i=0;i<outputValues.dimension(0);i++) {
212  for (int j=0;j<outputValues.dimension(1);j++) {
213  outputValues(i,j) = 0.0;
214  for (int k=0;k<this->getCardinality();k++) {
215  outputValues(i,j) += this->Vinv(k,i) * phisCur(k,j);
216  }
217  }
218  }
219  }
220  break;
221  case OPERATOR_GRAD:
222  case OPERATOR_D1:
223  case OPERATOR_D2:
224  case OPERATOR_D3:
225  case OPERATOR_D4:
226  case OPERATOR_D5:
227  case OPERATOR_D6:
228  case OPERATOR_D7:
229  case OPERATOR_D8:
230  case OPERATOR_D9:
231  case OPERATOR_D10:
232  {
233  const int dkcard =
234  (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,2): getDkCardinality(operatorType,2);
235 
236  FieldContainer<Scalar> phisCur( numBf , numPts , dkcard );
237  Phis.getValues( phisCur , inputPoints , operatorType );
238 
239  for (int i=0;i<outputValues.dimension(0);i++) {
240  for (int j=0;j<outputValues.dimension(1);j++) {
241  for (int k=0;k<outputValues.dimension(2);k++) {
242  outputValues(i,j,k) = 0.0;
243  for (int l=0;l<this->getCardinality();l++) {
244  outputValues(i,j,k) += this->Vinv(l,i) * phisCur(l,j,k);
245  }
246  }
247  }
248  }
249  }
250  break;
251  case OPERATOR_CURL: // only works in 2d. first component is -d/dy, second is d/dx
252  {
253  FieldContainer<Scalar> phisCur( numBf , numPts , getDkCardinality( OPERATOR_D1 , 2 ) );
254  Phis.getValues( phisCur , inputPoints , OPERATOR_D1 );
255 
256  for (int i=0;i<outputValues.dimension(0);i++) {
257  for (int j=0;j<outputValues.dimension(1);j++) {
258  outputValues(i,j,0) = 0.0;
259  outputValues(i,j,1) = 0.0;
260  for (int k=0;k<this->getCardinality();k++) {
261  outputValues(i,j,0) += this->Vinv(k,i) * phisCur(k,j,1);
262  }
263  for (int k=0;k<this->getCardinality();k++) {
264  outputValues(i,j,1) -= this->Vinv(k,i) * phisCur(k,j,0);
265  }
266  }
267  }
268  }
269  break;
270  default:
271  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
272  ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): Operator type not implemented");
273  break;
274  }
275  }
276  catch (std::invalid_argument &exception){
277  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
278  ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): Operator type not implemented");
279  }
280 
281  }
282 
283 
284 
285  template<class Scalar, class ArrayScalar>
286  void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
287  const ArrayScalar & inputPoints,
288  const ArrayScalar & cellVertices,
289  const EOperator operatorType) const {
290  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
291  ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): FEM Basis calling an FVD member function");
292  }
293 
294 
295 }// namespace Intrepid
296 #endif
297 
298 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
299 #ifdef __GNUC__
300 #warning "The Intrepid package is deprecated"
301 #endif
302 #endif
303 
Basis_HGRAD_TRI_Cn_FEM(const int n, const EPointType pointType)
Constructor.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell...