Intrepid
Intrepid_HCURL_QUAD_In_FEMDef.hpp
Go to the documentation of this file.
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
49namespace Intrepid {
50
51 template<class Scalar, class ArrayScalar>
53 const ArrayScalar & ptsClosed ,
54 const ArrayScalar & ptsOpen):
55 closedBasis_( order , ptsClosed ),
56 openBasis_( order-1 , ptsOpen ) ,
57 closedPts_( ptsClosed ),
58 openPts_( ptsOpen )
59 {
60 this -> basisDegree_ = order;
61 this -> basisCardinality_ = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality();
62 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
63 this -> basisType_ = BASIS_FEM_FIAT;
64 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
65 this -> basisTagsAreSet_ = false;
66
67 Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(2);
68 bases[0].resize(2); bases[1].resize(2);
69 bases[0][0] = rcp( &openBasis_ , false );
70 bases[0][1] = rcp( &closedBasis_ , false );
71 bases[1][0] = rcp( &closedBasis_ , false );
72 bases[1][1] = rcp( &openBasis_ , false );
73 this->setBases( bases );
74
75
76 }
77
78 template<class Scalar, class ArrayScalar>
80 closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ),
81 openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ) ,
82 closedPts_( order+1 , 1 ),
83 openPts_( order , 1 )
84 {
85 this -> basisDegree_ = order;
86 this -> basisCardinality_ = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality();
87 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
88 this -> basisType_ = BASIS_FEM_FIAT;
89 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
90 this -> basisTagsAreSet_ = false;
91
92 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ ,
93 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
94 order ,
95 0 ,
96 pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED );
97
98 if (pointType == POINTTYPE_SPECTRAL)
99 {
100 PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ ,
101 order - 1 );
102 }
103 else
104 {
105 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ ,
106 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
107 order + 1,
108 1 ,
109 POINTTYPE_EQUISPACED );
110 }
111
112 Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(2);
113 bases[0].resize(2); bases[1].resize(2);
114 bases[0][0] = rcp( &openBasis_ , false );
115 bases[0][1] = rcp( &closedBasis_ , false );
116 bases[1][0] = rcp( &closedBasis_ , false );
117 bases[1][1] = rcp( &openBasis_ , false );
118 this->setBases( bases );
119
120 }
121
122 template<class Scalar, class ArrayScalar>
124
125 // Basis-dependent intializations
126 int tagSize = 4; // size of DoF tag
127 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
128 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
129 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
130
131 std::vector<int> tags( tagSize * this->getCardinality() );
132
133 const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags();
134 const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags();
135
136 std::map<int,std::map<int,int> > total_dof_per_entity;
137 std::map<int,std::map<int,int> > current_dof_per_entity;
138
139 for (int i=0;i<4;i++) {
140 total_dof_per_entity[0][i] = 0;
141 current_dof_per_entity[0][i] = 0;
142 }
143 for (int i=0;i<4;i++) {
144 total_dof_per_entity[1][i] = 0;
145 current_dof_per_entity[1][i] = 0;
146 }
147 total_dof_per_entity[2][0] = 0;
148 current_dof_per_entity[2][0] = 0;
149
150 // tally dof on each facet. none on vertex
151 for (int i=0;i<4;i++) {
152 total_dof_per_entity[1][i] = openBasis_.getCardinality();
153 }
154
155 total_dof_per_entity[2][0] = this->getCardinality() - 4 * openBasis_.getCardinality();
156
157 int tagcur = 0;
158 // loop over the x-face tangent basis functions, which are (phi(x)psi(y),0)
159 // for psi in the closed basis and phi in the open
160 for (int j=0;j<closedBasis_.getCardinality();j++)
161 {
162 const int cdim = closedDofTags[j][0];
163 const int cent = closedDofTags[j][1];
164 for (int i=0;i<openBasis_.getCardinality();i++)
165 {
166 const int odim = openDofTags[i][0];
167 const int oent = openDofTags[i][1];
168 int dofdim;
169 int dofent;
170 ProductTopology::lineProduct2d(odim,oent,cdim,cent,dofdim,dofent);
171 tags[4*tagcur] = dofdim;
172 tags[4*tagcur+1] = dofent;
173 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
174 current_dof_per_entity[dofdim][dofent]++;
175 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
176 tagcur++;
177 }
178 }
179
180 // now we have to do it for the other basis functions, which are
181 // (psi(x)phi(y)) for psi in the closed basis and phi in the open
182 for (int j=0;j<openBasis_.getCardinality();j++)
183 {
184 const int odim = openDofTags[j][0];
185 const int oent = openDofTags[j][1];
186 for (int i=0;i<closedBasis_.getCardinality();i++)
187 {
188 const int cdim = closedDofTags[i][0];
189 const int cent = closedDofTags[i][1];
190 int dofdim;
191 int dofent;
192 ProductTopology::lineProduct2d(cdim,cent,odim,oent,dofdim,dofent);
193 tags[4*tagcur] = dofdim;
194 tags[4*tagcur+1] = dofent;
195 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
196 current_dof_per_entity[dofdim][dofent]++;
197 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
198 tagcur++;
199 }
200 }
201
202 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
203 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
204 this -> ordinalToTag_,
205 &(tags[0]),
206 this -> basisCardinality_,
207 tagSize,
208 posScDim,
209 posScOrd,
210 posDfOrd);
211 }
212
213
214 template<class Scalar, class ArrayScalar>
216 const ArrayScalar & inputPoints,
217 const EOperator operatorType) const {
218
219 // Verify arguments
220#ifdef HAVE_INTREPID_DEBUG
221 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
222 inputPoints,
223 operatorType,
224 this -> getBaseCellTopology(),
225 this -> getCardinality() );
226#endif
227
228 // Number of evaluation points = dim 0 of inputPoints
229 int dim0 = inputPoints.dimension(0);
230
231 // separate out points
232 FieldContainer<Scalar> xPoints(dim0,1);
233 FieldContainer<Scalar> yPoints(dim0,1);
234
235 for (int i=0;i<dim0;i++) {
236 xPoints(i,0) = inputPoints(i,0);
237 yPoints(i,0) = inputPoints(i,1);
238 }
239
240 switch (operatorType) {
241 case OPERATOR_VALUE:
242 {
243 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 );
244 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 );
245 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
246 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
247
248 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
249 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
250 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
251 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
252
253 int bfcur = 0;
254 // x component bfs are (open(x) closed(y),0)
255 for (int j=0;j<closedBasis_.getCardinality();j++)
256 {
257 for (int i=0;i<openBasis_.getCardinality();i++)
258 {
259 for (int l=0;l<dim0;l++)
260 {
261 outputValues(bfcur,l,0) = closedBasisValsYPts(j,l) * openBasisValsXPts(i,l);
262 outputValues(bfcur,l,1) = 0.0;
263 }
264 bfcur++;
265 }
266 }
267
268 // y component bfs are (0,closed(x) open(y))
269 for (int j=0;j<openBasis_.getCardinality();j++)
270 {
271 for (int i=0;i<closedBasis_.getCardinality();i++)
272 {
273 for (int l=0;l<dim0;l++)
274 {
275 outputValues(bfcur,l,0) = 0.0;
276 outputValues(bfcur,l,1) = openBasisValsYPts(j,l) * closedBasisValsXPts(i,l);
277 }
278 bfcur++;
279 }
280 }
281 }
282 break;
283 case OPERATOR_CURL:
284 {
285 FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 );
286 FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 );
287 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
288 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
289
290 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 );
291 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 );
292 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
293 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
294
295 int bfcur = 0;
296
297 // x component basis functions first
298 for (int j=0;j<closedBasis_.getCardinality();j++)
299 {
300 for (int i=0;i<openBasis_.getCardinality();i++)
301 {
302 for (int l=0;l<dim0;l++) {
303 outputValues(bfcur,l) = -closedBasisDerivsYPts(j,l,0) * openBasisValsXPts(i,l);
304 }
305 bfcur++;
306 }
307 }
308
309 // now y component basis functions
310 for (int j=0;j<openBasis_.getCardinality();j++)
311 {
312 for (int i=0;i<closedBasis_.getCardinality();i++)
313 {
314 for (int l=0;l<dim0;l++)
315 {
316 outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l);
317 }
318 bfcur++;
319 }
320 }
321 }
322 break;
323 case OPERATOR_DIV:
324 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
325 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): DIV is invalid operator for HCURL Basis Functions");
326 break;
327
328 case OPERATOR_GRAD:
329 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
330 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): GRAD is invalid operator for HCURL Basis Functions");
331 break;
332
333 case OPERATOR_D1:
334 case OPERATOR_D2:
335 case OPERATOR_D3:
336 case OPERATOR_D4:
337 case OPERATOR_D5:
338 case OPERATOR_D6:
339 case OPERATOR_D7:
340 case OPERATOR_D8:
341 case OPERATOR_D9:
342 case OPERATOR_D10:
343 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) ||
344 (operatorType == OPERATOR_D2) ||
345 (operatorType == OPERATOR_D3) ||
346 (operatorType == OPERATOR_D4) ||
347 (operatorType == OPERATOR_D5) ||
348 (operatorType == OPERATOR_D6) ||
349 (operatorType == OPERATOR_D7) ||
350 (operatorType == OPERATOR_D8) ||
351 (operatorType == OPERATOR_D9) ||
352 (operatorType == OPERATOR_D10) ),
353 std::invalid_argument,
354 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): Invalid operator type");
355 break;
356
357 default:
358 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
359 (operatorType != OPERATOR_GRAD) &&
360 (operatorType != OPERATOR_CURL) &&
361 (operatorType != OPERATOR_CURL) &&
362 (operatorType != OPERATOR_D1) &&
363 (operatorType != OPERATOR_D2) &&
364 (operatorType != OPERATOR_D3) &&
365 (operatorType != OPERATOR_D4) &&
366 (operatorType != OPERATOR_D5) &&
367 (operatorType != OPERATOR_D6) &&
368 (operatorType != OPERATOR_D7) &&
369 (operatorType != OPERATOR_D8) &&
370 (operatorType != OPERATOR_D9) &&
371 (operatorType != OPERATOR_D10) ),
372 std::invalid_argument,
373 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): Invalid operator type");
374 }
375 }
376
377
378
379 template<class Scalar, class ArrayScalar>
381 const ArrayScalar & inputPoints,
382 const ArrayScalar & cellVertices,
383 const EOperator operatorType) const {
384 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
385 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): FEM Basis calling an FVD member function");
386 }
387
388 template<class Scalar, class ArrayScalar>
390 {
391 // x-component basis functions
392 int cur = 0;
393
394 for (int j=0;j<closedPts_.dimension(0);j++)
395 {
396 for (int i=0;i<openPts_.dimension(0);i++)
397 {
398 DofCoords(cur,0) = openPts_(i,0);
399 DofCoords(cur,1) = closedPts_(j,0);
400 cur++;
401 }
402 }
403
404 // y-component basis functions
405 for (int j=0;j<openPts_.dimension(0);j++)
406 {
407 for (int i=0;i<closedPts_.dimension(0);i++)
408 {
409 DofCoords(cur,0) = closedPts_(i,0);
410 DofCoords(cur,1) = openPts_(j,0);
411 cur++;
412 }
413 }
414
415 return;
416 }
417
418
419
420}// 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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Quadrilateral cell.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
virtual void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
Basis_HCURL_QUAD_In_FEM(int order, const ArrayScalar &ptsClosed, const ArrayScalar &ptsOpen)
Constructor.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
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 lineProduct2d(const int dim0, const int entity0, const int dim1, const int entity1, int &resultdim, int &resultentity)