C++ CellJacobianBatch::jVals方法代码示例

void PeriodicSingleCellMesh1D::getJacobians(int cellDim, const Array<int>& cellLID,
    CellJacobianBatch& jBatch) const
{
  TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
    "cellDim=" << cellDim 
    << " is not in expected range [0, " << spatialDim()
    << "]");

  jBatch.resize(cellLID.size(), spatialDim(), cellDim);

  int nCells = cellLID.size();

  TEUCHOS_TEST_FOR_EXCEPT(nCells != 1);

  if (cellDim==0)
  {
    for (int i=0; i<nCells; i++)
    {
      double* detJ = jBatch.detJ(i);
      *detJ = 1.0;
    }
  }
  else
  {
    for (int i=0; i<nCells; i++)
    {
      double* J = jBatch.jVals(i);
      J[0] = fabs(xMin_-xMax_);
    }
  }
}

void CubicHermite::postApplyTransformationTriangle( const Mesh &mesh,
						    const Array<int> &cellLIDs , 
						    const CellJacobianBatch& JVol,
						    RCP<Array<double> >& A ) const
{
  Array<double> &Anoptr = *A;
  
  Array<double> cellVertH;
  getVertexH( mesh , cellLIDs , cellVertH );
  
  
  const int numRows = Anoptr.size() / 10 / JVol.numCells();
  
  for (int i=0;i<JVol.numCells();i++) 
    {
      const double *invJ = JVol.jVals(i);
      
      const int cell_start = i * numRows * 10;
      // handle columns 1 and 2
      for (int j=0;j<numRows;j++) 
	{
	  const double a = Anoptr[cell_start + numRows + j];
	  const double b = Anoptr[cell_start + 2*numRows + j];
	  Anoptr[cell_start + numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i];
	  Anoptr[cell_start + 2*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i];
	}
      
      // handle columns 4 and 5
      for (int j=0;j<numRows;j++) 
	{
	  const double a = Anoptr[cell_start + 4*numRows + j];
	  const double b = Anoptr[cell_start + 5*numRows + j];
	  Anoptr[cell_start + 4*numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i+1];
	  Anoptr[cell_start + 5*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i+1];
	}
      
      // handle columns 7 and 8
      for (int j=0;j<numRows;j++) 
	{
	  const double a = Anoptr[cell_start + 7*numRows + j];
	  const double b = Anoptr[cell_start + 8*numRows + j];
	  Anoptr[cell_start + 7*numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i+2];
	  Anoptr[cell_start + 8*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i+2];
	}
      
    }
  
  return;
}

void CubicHermite::preApplyTransformationTriangle( const Mesh &mesh, 
						   const Array<int> &cellLIDs,
						   const CellJacobianBatch& JVol,
						   RCP<Array<double> >& A) const
{
  Array<double> &Anoptr = *A;
  
  Array<double> cellVertH;
  getVertexH( mesh , cellLIDs , cellVertH );
  
  
  // this applies M from the left on each cell
  // A for each cell has 10 rows because it is Hermite
  // so, this gives the number of columns per cell
  const int numCols = Anoptr.size() / JVol.numCells() / 10;
  
  for (int i=0;i<JVol.numCells();i++) 
    {
      const int cell_start = i * numCols * 10;
      const double *invJ = JVol.jVals(i);
      
      for (int j=0;j<numCols;j++) 
	{
	  const int col_start = cell_start + 10 * j;
	  const double a1 = Anoptr[col_start + 1];
	  const double a2 = Anoptr[col_start + 2];
	  const double a4 = Anoptr[col_start + 4];
	  const double a5 = Anoptr[col_start + 5];
	  const double a7 = Anoptr[col_start + 7];
	  const double a8 = Anoptr[col_start + 8];
	  Anoptr[col_start+1] = (invJ[0]*a1 + invJ[1]*a2)/cellVertH[3*i];
	  Anoptr[col_start+2] = (invJ[2]*a1 + invJ[3]*a2)/cellVertH[3*i];
	  Anoptr[col_start+4] = (invJ[0]*a4 + invJ[1]*a5)/cellVertH[3*i+1];
	  Anoptr[col_start+5] = (invJ[2]*a4 + invJ[3]*a5)/cellVertH[3*i+1];
	  Anoptr[col_start+7] = (invJ[0]*a7 + invJ[1]*a8)/cellVertH[3*i+2];
	  Anoptr[col_start+8] = (invJ[2]*a7 + invJ[3]*a8)/cellVertH[3*i+2];
	}
    }
  
}

int main(int argc, char** argv)
{
  
  try
		{
      GlobalMPISession session(&argc, &argv);

      TimeMonitor t(totalTimer());

      int pMax = 2;
      int dim=2;

      verbosity<RefIntegral>() = 0;

      CellType cellType = TriangleCell;

      Point a = Point(0.0, 0.0);
      Point b = Point(1.0, 0.0);
      Point c = Point(0.0, 1.0);

//       Point a = Point(1.0, 1.0);
//       Point b = Point(1.2, 1.6);
//       Point c = Point(0.8, 1.3);
      CellJacobianBatch JBatch;
      JBatch.resize(1, dim, dim);
      double* J = JBatch.jVals(0);
      J[0] = b[0] - a[0];
      J[1] = c[0] - a[0];
      J[2] = b[1] - a[1];
      J[3] = c[1] - a[1];

      QuadratureFamily q4 = new GaussianQuadrature(4);
          
      for (int p=1; p<=pMax; p++)
        {
          std::cerr << std::endl << "---------- p = " << p << " --------------" << std::endl;
          Tabs tab;
          BasisFamily P = new Lagrange(p);
      
          RCP<Array<double> > A = rcp(new Array<double>());
          RCP<Array<double> > Bxx = rcp(new Array<double>());
          RCP<Array<double> > Byy = rcp(new Array<double>());

          Array<double> constCoeff = tuple(1.0, 1.0);

          Array<int> alpha = tuple(0,1);
          Array<int> beta = tuple(0,1);
          ParametrizedCurve curve = new DummyParametrizedCurve();
           MeshType meshType = new BasicSimplicialMeshType();
           MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, 10, meshType);
           Mesh mesh = mesher.getMesh();
           QuadratureFamily quad_1 = new GaussianQuadrature(2);
           RCP<Array<int> > cellLIDs;

          RefIntegral ref(dim, dim, cellType, P, alpha, 1, P, beta,quad_1, 1, isInternalBdry , curve , mesh);
          QuadratureIntegral qxx(dim, dim, cellType, P, tuple(0), 1, 
                                 P, tuple(0), 1, q4, curve , mesh,isInternalBdry);
          QuadratureIntegral qyy(dim, dim, cellType, P, tuple(1), 1, 
                                 P, tuple(1), 1, q4, curve , mesh,isInternalBdry);

          int nq = qxx.nQuad();
          Array<double> varCoeff(nq, 1.0);

          std::cerr << "============================= Doing reference integration =========== " << std::endl;

          ref.transformTwoForm(JBatch, constCoeff, cellLIDs, A);
          std::cerr << "============================= Doing quad integration xx =========== " << std::endl;
          qxx.transformTwoForm(JBatch, &(varCoeff[0]), cellLIDs , Bxx);
          std::cerr << "============================= Doing quad integration yy =========== " << std::endl;
          qyy.transformTwoForm(JBatch, &(varCoeff[0]), cellLIDs , Byy);

          std::cerr << "============================= Done integration =========== " << std::endl;
          std::cerr << tab << "transformed reference element" << std::endl;
          std::cerr << tab << "{";
          for (int r=0; r<ref.nNodesTest(); r++)
            {
              if (r!=0) std::cerr << ", ";
              std::cerr << "{";
              for (int c=0; c<ref.nNodesUnk(); c++)
                {
                  if (c!=0) std::cerr << ", ";
                  std::cerr << (*A)[r + ref.nNodesTest()*c];
                }
              std::cerr << "}";
            }
          std::cerr << "}" << std::endl;

          std::cerr << tab << "transformed Q_xx" << std::endl;
          std::cerr << tab << "{";
          for (int r=0; r<qxx.nNodesTest(); r++)
            {
              if (r!=0) std::cerr << ", ";
              std::cerr << "{";
              for (int c=0; c<qxx.nNodesUnk(); c++)
                {
                  if (c!=0) std::cerr << ", ";
                  int i = r + qxx.nNodesTest()*c;
                  double lapl = (*Bxx)[i];
                  std::cerr << lapl;
                }
//.........这里部分代码省略.........

void PeanoMesh3D::getJacobians(int cellDim, const Array<int>& cellLID,
                          CellJacobianBatch& jBatch) const
{
	  //printf("cellDim:%d  _uniqueResolution:%f ",cellDim, _uniqueResolution);
	  SUNDANCE_VERB_HIGH("getJacobians()");
	  TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
	    "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
	  int nCells = cellLID.size();
 	  int tmp_index , tmp;
 	  int tmp_index1 , tmp_index2;
 	  Point pnt(0.0,0.0,0.0);
 	  Point pnt1(0.0,0.0,0.0);
 	  Point pnt2(0.0,0.0,0.0);
	  jBatch.resize(cellLID.size(), spatialDim(), cellDim);
	  if (cellDim < spatialDim()) // they need the Jacobian of a lower dinemsional element
	  {
		  //printf("PeanoMesh3D::getJacobians() cellDim < spatialDim() \n");
		   for (int i=0; i<nCells; i++)
		    {
		      //printf("PeanoMesh3D::getJacobian() cellDim < spatialDim() cellDim:%d , ret:%f \n",cellDim , _uniqueResolution);
		      double* detJ = jBatch.detJ(i);
		      switch(cellDim)
		      {
		        case 0: *detJ = 1.0;
		          break;
		        case 1:
			    	 tmp_index = this->facetLID(cellDim,  cellLID[i] , 0 , 0 , tmp );
			    	 tmp_index1= this->facetLID(cellDim,  cellLID[i] , 0 , 1 , tmp );
			    	 pnt = nodePosition(tmp_index);
			    	 pnt1 = nodePosition(tmp_index1);
			    	 pnt = pnt1 - pnt;
		             *detJ = sqrt(pnt * pnt); // the length of the edge
		        break;
		        case 2:{
			    	 tmp_index = this->facetLID(cellDim,  cellLID[i] , 0 , 0 , tmp );
			    	 tmp_index1= this->facetLID(cellDim,  cellLID[i] , 0 , 1 , tmp );
			    	 tmp_index2= this->facetLID(cellDim,  cellLID[i] , 0 , 2 , tmp );
			    	 pnt = nodePosition(tmp_index);
			    	 pnt1 = nodePosition(tmp_index1);
			    	 pnt2 = nodePosition(tmp_index2);
			    	 Point directedArea = cross( pnt1 - pnt , pnt2 - pnt );
		             *detJ = sqrt(directedArea * directedArea); // the are of the face
		        break;}
		        default:
		          TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
		            "cellDim=" << cellDim << " in PeanoMesh3D::getJacobians()");
		      }
		    }
	  }else{ // they request the complete Jacoby matrix for this bunch of elements
		    //Array<double> J(cellDim*cellDim);
		    SUNDANCE_VERB_HIGH("cellDim == spatialDim()");
		    for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
		    {
			  //printf("PeanoMesh3D::getJacobian() cellDim == spatialDim() cellDim:%d , ret:%f \n",cellDim , _uniqueResolution);
		      double* J = jBatch.jVals(i);
		      switch(cellDim)
		      {
		        case 3:
		          J[0] = _peanoMesh->returnResolution(0);
		          J[1] = 0.0; J[2] = 0.0; J[3] = 0.0;
		          J[4] = _peanoMesh->returnResolution(1);
		          J[5] = 0.0; J[6] = 0.0; J[7] = 0.0;
		          J[8] = _peanoMesh->returnResolution(2); // the Jacobi of the tet
		        break;
		        default:
		          TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
		            "cellDim=" << cellDim
		            << " in PeanoMesh3D::getJacobians()");
		      }
		    }
	  }
}

int main(int argc, char** argv)
{
  
  try
  {
    GlobalMPISession session(&argc, &argv);

    TimeMonitor t(totalTimer());

    int pMax = 1;
    int dim=2;

    CellType cellType = TriangleCell;

    Point a = Point(0.0, 0.0);
    Point b = Point(1.0, 0.0);
    Point c = Point(0.0, 1.0);
    CellJacobianBatch JBatch;
    JBatch.resize(1, 2, 2);
    double* J = JBatch.jVals(0);
    J[0] = b[0] - a[0];
    J[1] = c[0] - a[0];
    J[2] = b[1] - a[1];
    J[3] = c[1] - a[1];


    bool isInternalBdry=false;

    /* ------ evaluate Lagrange and FIAT-Lagrange at the vertices */
    Array<Point> verts = tuple(a,b,c);
    BasisFamily lagrange = new Lagrange(1);
    BasisFamily fiatLagrange = new Lagrange(1);
      
    MultiIndex d0(0,0,0);
    MultiIndex dx(1,0,0);
    MultiIndex dy(0,1,0);

    Array<Array<Array<double> > > result;

    Array<int> dummy;

    std::cerr << "------ Evaluating bases at vertices ----------" << std::endl
         << std::endl;

    std::cerr << "Evaluating phi(vert) with FIAT-Lagrange" << std::endl;
    fiatLagrange.ptr()->refEval(cellType, verts, d0, result);
    std::cerr << "results = " << result << std::endl << std::endl;

    std::cerr << "Evaluating phi(vert) with Lagrange" << std::endl;
    lagrange.ptr()->refEval(cellType, verts, d0, result);
    std::cerr << "results = " << result << std::endl << std::endl;

    std::cerr << std::endl ;

    std::cerr << "Evaluating Dx*phi(vert) with FIAT-Lagrange" << std::endl;
    fiatLagrange.ptr()->refEval(cellType, verts, dx, result);
    std::cerr << "results = " << result << std::endl << std::endl;

    std::cerr << "Evaluating Dx*phi(vert) with Lagrange" << std::endl;
    lagrange.ptr()->refEval(cellType, verts, dx, result);
    std::cerr << "results = " << result << std::endl << std::endl;

    std::cerr << std::endl ;
      
    std::cerr << "Evaluating Dy*phi(vert) with FIAT-Lagrange" << std::endl;
    fiatLagrange.ptr()->refEval(cellType, verts, dy, result);
    std::cerr << "results = " << result << std::endl << std::endl;

    std::cerr << "Evaluating Dy*phi(vert) with Lagrange" << std::endl;
    lagrange.ptr()->refEval(cellType, verts, dy, result);
    std::cerr << "results = " << result << std::endl << std::endl;

      

    /* --------- evaluate integrals over elements ----------- */
      
    RCP<Array<double> > A = rcp(new Array<double>());
          
    QuadratureFamily quad = new GaussianQuadrature(4);
    Array<double> quadWeights;
    Array<Point> quadPts;
    quad.getPoints(cellType, quadPts, quadWeights);
    int nQuad = quadPts.size();

    Array<double> coeff(nQuad);
    for (int i=0; i<nQuad; i++) 
    {
      double s = quadPts[i][0];
      double t = quadPts[i][1];
      double x = a[0] + J[0]*s + J[1]*t;
      double y = a[1] + J[2]*s + J[3]*t;
      coeff[i] = x*y;
    }
    const double* const f = &(coeff[0]);

    std::cerr << std::endl << std::endl 
         << "---------------- One-forms --------------------" 
         << std::endl << std::endl;
    for (int p=1; p<=pMax; p++)
    {
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  int stat = 0;
  int verb=1;
  try
  {
    GlobalMPISession session(&argc, &argv);

    TimeMonitor t(totalTimer());

    int pMax = 2;
    int dim=2;

    bool isInternalBdry = false;

    Utils::setChopVal(1.0e-14);

    CellType cellType = TriangleCell;

    //       Point a = Point(1.0, 1.0);
    //       Point b = Point(1.2, 1.6);
    //       Point c = Point(0.8, 1.3);

    Point a = Point(0.0, 0.0);
    Point b = Point(1.0, 0.0);
    Point c = Point(0.0, 1.0);

    Point d = Point(0.1, 0.1);
    Point e = Point(1.0, 0.0);
    Point f = Point(0.0, 1.0);

    int nCells = 2;

    CellJacobianBatch JBatch;
    JBatch.resize(nCells, 2, 2);
    double* J = JBatch.jVals(0);
    J[0] = b[0] - a[0];
    J[1] = c[0] - a[0];
    J[2] = b[1] - a[1];
    J[3] = c[1] - a[1];

    J[4] = e[0] - d[0];
    J[5] = f[0] - d[0];
    J[6] = e[1] - d[1];
    J[7] = f[1] - d[1];


      
    Array<int> dummy;
    double coeff = 1.0;
    RCP<Array<double> > A = rcp(new Array<double>());
    RCP<Array<double> > B = rcp(new Array<double>());

    QuadratureFamily q4 = new GaussianQuadrature(4);

    int nErrors = 0;

    std::cerr << std::endl << std::endl 
         << "---------------- One-forms --------------------" 
         << std::endl << std::endl;
    for (int p=0; p<=pMax; p++)
    {
      BasisFamily P = new Lagrange(p);
      for (int dp=0; dp<=1; dp++)
      {
        if (dp > p) continue;

        int numTestDir = 1;
        if (dp==1) numTestDir = dim;
        for (int t=0; t<numTestDir; t++)
        {
          int alpha = t;
          Tabs tab;

          ParametrizedCurve curve = new DummyParametrizedCurve();
          MeshType meshType = new BasicSimplicialMeshType();
          MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, 10, meshType);
          Mesh mesh = mesher.getMesh();
          RCP<Array<int> > cellLIDs;

          RefIntegral ref(dim, cellType, dim, cellType, P, alpha, dp, q4 , isInternalBdry, curve, mesh ,verb);
          A->resize(JBatch.numCells() * ref.nNodes());
          for (int ai=0; ai<A->size(); ai++) (*A)[ai]=0.0;
          ref.transformOneForm(JBatch, JBatch, dummy, cellLIDs , coeff, A);
          std::cerr << tab << "transformed reference element" << std::endl;
          if (dp>0) std::cerr << tab << "test diff direction=" << t << std::endl;
          for (int cell=0; cell<nCells; cell++)
          {
            std::cerr << tab << "{";
            for (int r=0; r<ref.nNodesTest(); r++)
            {
              if (r!=0) std::cerr << ", ";
              std::cerr << Utils::chop((*A)[cell*ref.nNodesTest()+r]);
            }
            std::cerr << "}" << std::endl;
          }
          QuadratureIntegral quad(dim, cellType, dim, cellType, P, alpha, dp, q4, isInternalBdry, curve, mesh, verb);
          Array<double> quadCoeff(2*quad.nQuad(), 1.0);
          B->resize(JBatch.numCells() * quad.nNodes());
          for (int ai=0; ai<B->size(); ai++) (*B)[ai]=0.0;
//.........这里部分代码省略.........