Albany: a Trilinos-based PDE code: SurfaceBasis

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 #include <Intrepid_MiniTensor.h>
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009 
00010 #include "Sacado_MathFunctions.hpp"
00011 
00012 namespace LCM {
00013 
00014 //----------------------------------------------------------------------
00015   template<typename EvalT, typename Traits>
00016   SurfaceBasis<EvalT, Traits>::SurfaceBasis(const Teuchos::ParameterList& p,
00017                                             const Teuchos::RCP<Albany::Layouts>& dl) :
00018       needCurrentBasis(false),
00019       referenceCoords(p.get<std::string>("Reference Coordinates Name"), dl->vertices_vector),
00020       cubature       (p.get<Teuchos::RCP<Intrepid::Cubature<RealType> > >("Cubature")),
00021       intrepidBasis  (p.get<Teuchos::RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > > >("Intrepid Basis")),
00022       refBasis       (p.get<std::string>("Reference Basis Name"), dl->qp_tensor),
00023       refArea        (p.get<std::string>("Reference Area Name"), dl->qp_scalar),
00024       refDualBasis   (p.get<std::string>("Reference Dual Basis Name"), dl->qp_tensor),
00025       refNormal      (p.get<std::string>("Reference Normal Name"), dl->qp_vector)
00026   {
00027     this->addDependentField(referenceCoords);
00028     this->addEvaluatedField(refBasis);
00029     this->addEvaluatedField(refArea);
00030     this->addEvaluatedField(refDualBasis);
00031     this->addEvaluatedField(refNormal);
00032 
00033     // if current coordinates are being passed in, compute and return the current basis
00034     // needed for the localization element, but not uncoupled transport
00035     if (p.isType<std::string>("Current Coordinates Name")) {
00036       needCurrentBasis = true;
00037 
00038       // grab the current coords
00039       PHX::MDField<ScalarT, Cell, Vertex, Dim> tmp(p.get<std::string>("Current Coordinates Name"), dl->node_vector);
00040       currentCoords = tmp;
00041 
00042       // set up the current basis
00043       PHX::MDField<ScalarT, Cell, QuadPoint, Dim, Dim> tmp2(p.get<std::string>("Current Basis Name"), dl->qp_tensor);
00044       currentBasis = tmp2;
00045 
00046       this->addDependentField(currentCoords);
00047       this->addEvaluatedField(currentBasis);
00048     }
00049 
00050     // Get Dimensions
00051     std::vector<PHX::DataLayout::size_type> dims;
00052     dl->node_vector->dimensions(dims);
00053 
00054     int containerSize = dims[0];
00055     numNodes = dims[1];
00056     numPlaneNodes = numNodes / 2;
00057 
00058     numQPs = cubature->getNumPoints();
00059     numPlaneDims = cubature->getDimension();
00060     numDims = numPlaneDims + 1;
00061 
00062 #ifdef ALBANY_VERBOSE
00063     std::cout << "in Surface Basis" << std::endl;
00064     std::cout << " numPlaneNodes: " << numPlaneNodes << std::endl;
00065     std::cout << " numPlaneDims: " << numPlaneDims << std::endl;
00066     std::cout << " numQPs: " << numQPs << std::endl;
00067     std::cout << " cubature->getNumPoints(): " << cubature->getNumPoints() << std::endl;
00068     std::cout << " cubature->getDimension(): " << cubature->getDimension() << std::endl;
00069 #endif
00070 
00071     // Allocate Temporary FieldContainers
00072     refValues.resize(numPlaneNodes, numQPs);
00073     refGrads.resize(numPlaneNodes, numQPs, numPlaneDims);
00074     refPoints.resize(numQPs, numPlaneDims);
00075     refWeights.resize(numQPs);
00076 
00077     // temp space for midplane coords
00078     refMidplaneCoords.resize(containerSize, numPlaneNodes, numDims);
00079     currentMidplaneCoords.resize(containerSize, numPlaneNodes, numDims);
00080 
00081     // Pre-Calculate reference element quantitites
00082     cubature->getCubature(refPoints, refWeights);
00083     intrepidBasis->getValues(refValues, refPoints, Intrepid::OPERATOR_VALUE);
00084     intrepidBasis->getValues(refGrads, refPoints, Intrepid::OPERATOR_GRAD);
00085 
00086     this->setName("SurfaceBasis" + PHX::TypeString<EvalT>::value);
00087   }
00088 
00089   //----------------------------------------------------------------------
00090   template<typename EvalT, typename Traits>
00091   void SurfaceBasis<EvalT, Traits>::postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager<Traits>& fm)
00092   {
00093     this->utils.setFieldData(referenceCoords, fm);
00094     this->utils.setFieldData(refArea, fm);
00095     this->utils.setFieldData(refDualBasis, fm);
00096     this->utils.setFieldData(refNormal, fm);
00097     this->utils.setFieldData(refBasis, fm);
00098     if (needCurrentBasis) {
00099       this->utils.setFieldData(currentCoords, fm);
00100       this->utils.setFieldData(currentBasis, fm);
00101     }
00102   }
00103 
00104   //----------------------------------------------------------------------
00105   template<typename EvalT, typename Traits>
00106   void SurfaceBasis<EvalT, Traits>::evaluateFields(typename Traits::EvalData workset)
00107   {
00108     for (std::size_t cell(0); cell < workset.numCells; ++cell) {
00109       // for the reference geometry
00110       // compute the mid-plane coordinates
00111       computeReferenceMidplaneCoords(referenceCoords, refMidplaneCoords);
00112 
00113       // compute basis vectors
00114       computeReferenceBaseVectors(refMidplaneCoords, refBasis);
00115 
00116       // compute the dual
00117       computeDualBaseVectors(refMidplaneCoords, refBasis, refNormal, refDualBasis);
00118 
00119       // compute the Jacobian
00120       computeJacobian(refBasis, refDualBasis, refArea);
00121 
00122       if (needCurrentBasis) {
00123         // for the current configuration
00124         // compute the mid-plane coordinates
00125         computeCurrentMidplaneCoords(currentCoords, currentMidplaneCoords);
00126 
00127         // compute base vectors
00128         computeCurrentBaseVectors(currentMidplaneCoords, currentBasis);
00129       }
00130     }
00131   }
00132   //----------------------------------------------------------------------
00133   template<typename EvalT, typename Traits>
00134   void SurfaceBasis<EvalT, Traits>::computeReferenceMidplaneCoords(PHX::MDField<MeshScalarT, Cell, Vertex, Dim> coords,
00135                                                                    MFC & midplaneCoords)
00136   {
00137     for (int cell(0); cell < midplaneCoords.dimension(0); ++cell) {
00138       // compute the mid-plane coordinates
00139       for (int node(0); node < numPlaneNodes; ++node) {
00140         int topNode = node + numPlaneNodes;
00141         for (int dim(0); dim < numDims; ++dim) {
00142           midplaneCoords(cell, node, dim) = 0.5 * (coords(cell, node, dim) + coords(cell, topNode, dim));
00143         }
00144       }
00145     }
00146   }
00147   //----------------------------------------------------------------------
00148   template<typename EvalT, typename Traits>
00149   void SurfaceBasis<EvalT, Traits>::computeCurrentMidplaneCoords(PHX::MDField<ScalarT, Cell, Vertex, Dim> coords,
00150                                                                  SFC & midplaneCoords)
00151   {
00152     for (int cell(0); cell < midplaneCoords.dimension(0); ++cell) {
00153       // compute the mid-plane coordinates
00154       for (int node(0); node < numPlaneNodes; ++node) {
00155         int topNode = node + numPlaneNodes;
00156         for (int dim(0); dim < numDims; ++dim) {
00157           midplaneCoords(cell, node, dim) = 0.5 * (coords(cell, node, dim) + coords(cell, topNode, dim));
00158         }
00159       }
00160     }
00161   }
00162   //----------------------------------------------------------------------
00163   template<typename EvalT, typename Traits>
00164   void SurfaceBasis<EvalT, Traits>::
00165   computeReferenceBaseVectors(const MFC & midplaneCoords,
00166                               PHX::MDField<MeshScalarT, Cell, QuadPoint, Dim, Dim> basis)
00167   {
00168     for (int cell(0); cell < midplaneCoords.dimension(0); ++cell) {
00169       // get the midplane coordinates
00170       std::vector<Intrepid::Vector<MeshScalarT> > midplaneNodes(numPlaneNodes);
00171       for (std::size_t node(0); node < numPlaneNodes; ++node)
00172         midplaneNodes[node] = Intrepid::Vector<MeshScalarT>(3, &midplaneCoords(cell, node, 0));
00173 
00174       Intrepid::Vector<MeshScalarT> g_0(0, 0, 0), g_1(0, 0, 0), g_2(0, 0, 0);
00175       //compute the base vectors
00176       for (std::size_t pt(0); pt < numQPs; ++pt) {
00177         g_0.clear();
00178         g_1.clear();
00179         g_2.clear();
00180         for (std::size_t node(0); node < numPlaneNodes; ++node) {
00181           g_0 += refGrads(node, pt, 0) * midplaneNodes[node];
00182           g_1 += refGrads(node, pt, 1) * midplaneNodes[node];
00183         }
00184         g_2 = cross(g_0, g_1) / norm(cross(g_0, g_1));
00185 
00186         basis(cell, pt, 0, 0) = g_0(0);
00187         basis(cell, pt, 0, 1) = g_0(1);
00188         basis(cell, pt, 0, 2) = g_0(2);
00189         basis(cell, pt, 1, 0) = g_1(0);
00190         basis(cell, pt, 1, 1) = g_1(1);
00191         basis(cell, pt, 1, 2) = g_1(2);
00192         basis(cell, pt, 2, 0) = g_2(0);
00193         basis(cell, pt, 2, 1) = g_2(1);
00194         basis(cell, pt, 2, 2) = g_2(2);
00195       }
00196     }
00197   }
00198   //----------------------------------------------------------------------
00199   template<typename EvalT, typename Traits>
00200   void SurfaceBasis<EvalT, Traits>::
00201   computeCurrentBaseVectors(const SFC & midplaneCoords,
00202                             PHX::MDField<ScalarT, Cell, QuadPoint, Dim, Dim> basis)
00203   {
00204     for (int cell(0); cell < midplaneCoords.dimension(0); ++cell) {
00205       // get the midplane coordinates
00206       std::vector<Intrepid::Vector<ScalarT> > midplaneNodes(numPlaneNodes);
00207       for (std::size_t node(0); node < numPlaneNodes; ++node)
00208         midplaneNodes[node] = Intrepid::Vector<ScalarT>(3, &midplaneCoords(cell, node, 0));
00209 
00210       Intrepid::Vector<ScalarT> g_0(0, 0, 0), g_1(0, 0, 0), g_2(0, 0, 0);
00211       //compute the base vectors
00212       for (std::size_t pt(0); pt < numQPs; ++pt) {
00213         g_0.clear();
00214         g_1.clear();
00215         g_2.clear();
00216         for (std::size_t node(0); node < numPlaneNodes; ++node) {
00217           g_0 += refGrads(node, pt, 0) * midplaneNodes[node];
00218           g_1 += refGrads(node, pt, 1) * midplaneNodes[node];
00219         }
00220         g_2 = cross(g_0, g_1) / norm(cross(g_0, g_1));
00221 
00222         basis(cell, pt, 0, 0) = g_0(0);
00223         basis(cell, pt, 0, 1) = g_0(1);
00224         basis(cell, pt, 0, 2) = g_0(2);
00225         basis(cell, pt, 1, 0) = g_1(0);
00226         basis(cell, pt, 1, 1) = g_1(1);
00227         basis(cell, pt, 1, 2) = g_1(2);
00228         basis(cell, pt, 2, 0) = g_2(0);
00229         basis(cell, pt, 2, 1) = g_2(1);
00230         basis(cell, pt, 2, 2) = g_2(2);
00231       }
00232     }
00233   }
00234   //----------------------------------------------------------------------
00235   template<typename EvalT, typename Traits>
00236   void SurfaceBasis<EvalT, Traits>::computeDualBaseVectors(
00237       const MFC & midplaneCoords,
00238       const PHX::MDField<MeshScalarT, Cell, QuadPoint, Dim, Dim> basis,
00239       PHX::MDField<MeshScalarT, Cell, QuadPoint, Dim> normal,
00240       PHX::MDField<MeshScalarT, Cell, QuadPoint, Dim, Dim> dualBasis)
00241   {
00242     std::size_t worksetSize = midplaneCoords.dimension(0);
00243 
00244     Intrepid::Vector<MeshScalarT> g_0(0, 0, 0), g_1(0, 0, 0), g_2(0, 0, 0), g0(0, 0, 0),
00245         g1(0, 0, 0), g2(0, 0, 0);
00246 
00247     for (std::size_t cell(0); cell < worksetSize; ++cell) {
00248       for (std::size_t pt(0); pt < numQPs; ++pt) {
00249         g_0 = Intrepid::Vector<MeshScalarT>(3, &basis(cell, pt, 0, 0));
00250         g_1 = Intrepid::Vector<MeshScalarT>(3, &basis(cell, pt, 1, 0));
00251         g_2 = Intrepid::Vector<MeshScalarT>(3, &basis(cell, pt, 2, 0));
00252 
00253         normal(cell, pt, 0) = g_2(0);
00254         normal(cell, pt, 1) = g_2(1);
00255         normal(cell, pt, 2) = g_2(2);
00256 
00257         g0 = cross(g_1, g_2) / dot(g_0, cross(g_1, g_2));
00258         g1 = cross(g_0, g_2) / dot(g_1, cross(g_0, g_2));
00259         g2 = cross(g_0, g_1) / dot(g_2, cross(g_0, g_1));
00260 
00261         dualBasis(cell, pt, 0, 0) = g0(0);
00262         dualBasis(cell, pt, 0, 1) = g0(1);
00263         dualBasis(cell, pt, 0, 2) = g0(2);
00264         dualBasis(cell, pt, 1, 0) = g1(0);
00265         dualBasis(cell, pt, 1, 1) = g1(1);
00266         dualBasis(cell, pt, 1, 2) = g1(2);
00267         dualBasis(cell, pt, 2, 0) = g2(0);
00268         dualBasis(cell, pt, 2, 1) = g2(1);
00269         dualBasis(cell, pt, 2, 2) = g2(2);
00270       }
00271     }
00272   }
00273   //----------------------------------------------------------------------
00274   template<typename EvalT, typename Traits>
00275   void SurfaceBasis<EvalT, Traits>::computeJacobian(
00276       const PHX::MDField<MeshScalarT, Cell, QuadPoint, Dim, Dim> basis,
00277       const PHX::MDField<MeshScalarT, Cell, QuadPoint, Dim, Dim> dualBasis,
00278       PHX::MDField<MeshScalarT, Cell, QuadPoint> area)
00279   {
00280     const std::size_t worksetSize = basis.dimension(0);
00281 
00282     for (std::size_t cell(0); cell < worksetSize; ++cell) {
00283       for (std::size_t pt(0); pt < numQPs; ++pt) {
00284         Intrepid::Tensor<MeshScalarT> dPhiInv(3, &dualBasis(cell, pt, 0, 0));
00285         Intrepid::Tensor<MeshScalarT> dPhi(3, &basis(cell, pt, 0, 0));
00286         Intrepid::Vector<MeshScalarT> G_2(3, &basis(cell, pt, 2, 0));
00287 
00288         MeshScalarT j0 = Intrepid::det(dPhi);
00289         MeshScalarT jacobian = j0 *
00290           std::sqrt( Intrepid::dot(Intrepid::dot(G_2, Intrepid::transpose(dPhiInv) * dPhiInv), G_2));
00291         area(cell, pt) = jacobian * refWeights(pt);
00292       }
00293     }
00294 
00295   }
00296 //----------------------------------------------------------------------
00297 }//namespace LCM
SurfaceBasis_Def.hpp
