C++ BasisCachePtr::cellIDs方法代码示例

  virtual void values(FieldContainer<double> &values, BasisCachePtr basisCache)
  {
    // not the most efficient implementation
    _fxn->values(values,basisCache);
    vector<GlobalIndexType> contextCellIDs = basisCache->cellIDs();
    int cellIndex=0; // keep track of index into values

    int entryCount = values.size();
    int numCells = values.dimension(0);
    int numEntriesPerCell = entryCount / numCells;

    for (vector<GlobalIndexType>::iterator cellIt = contextCellIDs.begin(); cellIt != contextCellIDs.end(); cellIt++)
    {
      GlobalIndexType cellID = *cellIt;
      if (_cellIDs.find(cellID) == _cellIDs.end())
      {
        // clear out the associated entries
        for (int j=0; j<numEntriesPerCell; j++)
        {
          values[cellIndex*numEntriesPerCell + j] = 0;
        }
      }
      cellIndex++;
    }
  }

bool checkLTSumConsistency(LinearTermPtr a, LinearTermPtr b, DofOrderingPtr dofOrdering, BasisCachePtr basisCache)
{
  double tol = 1e-14;

  int numCells = basisCache->cellIDs().size();
  int numDofs = dofOrdering->totalDofs();
  bool forceBoundaryTerm = false;
  FieldContainer<double> aValues(numCells,numDofs), bValues(numCells,numDofs), sumValues(numCells,numDofs);
  a->integrate(aValues,dofOrdering,basisCache,forceBoundaryTerm);
  b->integrate(bValues,dofOrdering,basisCache,forceBoundaryTerm);
  (a+b)->integrate(sumValues, dofOrdering, basisCache, forceBoundaryTerm);

  int size = aValues.size();

  for (int i=0; i<size; i++)
  {
    double expectedValue = aValues[i] + bValues[i];
    double diff = abs( expectedValue - sumValues[i] );
    if (diff > tol)
    {
      return false;
    }
  }
  return true;
}

 void values(FieldContainer<double> &values, BasisCachePtr sliceBasisCache) {
   vector<GlobalIndexType> sliceCellIDs = sliceBasisCache->cellIDs();
   
   Teuchos::Array<int> dim;
   values.dimensions(dim);
   dim[0] = 1; // one cell
   Teuchos::Array<int> offset(dim.size());
   
   for (int cellOrdinal = 0; cellOrdinal < sliceCellIDs.size(); cellOrdinal++) {
     offset[0] = cellOrdinal;
     int enumeration = values.getEnumeration(offset);
     FieldContainer<double>valuesForCell(dim,&values[enumeration]);
     GlobalIndexType sliceCellID = sliceCellIDs[cellOrdinal];
     int numPoints = sliceBasisCache->getPhysicalCubaturePoints().dimension(1);
     int spaceDim = sliceBasisCache->getPhysicalCubaturePoints().dimension(2);
     FieldContainer<double> spaceTimePhysicalPoints(1,numPoints,spaceDim+1);
     for (int ptOrdinal=0; ptOrdinal<numPoints; ptOrdinal++) {
       for (int d=0; d<spaceDim; d++) {
         spaceTimePhysicalPoints(0,ptOrdinal,d) = sliceBasisCache->getPhysicalCubaturePoints()(cellOrdinal,ptOrdinal,d);
       }
       spaceTimePhysicalPoints(0,ptOrdinal,spaceDim) = _t;
     }
     
     GlobalIndexType cellID = _cellIDMap[sliceCellID];
     BasisCachePtr spaceTimeBasisCache = BasisCache::basisCacheForCell(_spaceTimeMesh, cellID);
     
     FieldContainer<double> spaceTimeRefPoints(1,numPoints,spaceDim+1);
     CamelliaCellTools::mapToReferenceFrame(spaceTimeRefPoints, spaceTimePhysicalPoints, _spaceTimeMesh, cellID);
     spaceTimeRefPoints.resize(numPoints,spaceDim+1);
     spaceTimeBasisCache->setRefCellPoints(spaceTimeRefPoints);
     _spaceTimeFunction->values(valuesForCell, spaceTimeBasisCache);
   }
 }

  void values(FieldContainer<double> &values, BasisCachePtr basisCache)
  {
    int numCells = values.dimension(0);
    int numPoints = values.dimension(1);

    MeshPtr mesh = basisCache->mesh();
    vector<int> cellIDs = basisCache->cellIDs();

    double tol=1e-14;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      double h = 1.0;
      if (_spatialCoord==0)
      {
        h = mesh->getCellXSize(cellIDs[cellIndex]);
      }
      else if (_spatialCoord==1)
      {
        h = mesh->getCellYSize(cellIDs[cellIndex]);
      }
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        values(cellIndex,ptIndex) = sqrt(h);
      }
    }
  }

bool FunctionTests::testAdaptiveIntegrate()
{
  bool success = true;

  // we must create our own basisCache here because _basisCache
  // has had its ref cell points set, which basically means it's
  // opted out of having any help with integration.
  BasisCachePtr basisCache = Teuchos::rcp( new BasisCache( _elemType, _spectralConfusionMesh ) );
  vector<GlobalIndexType> cellIDs;
  cellIDs.push_back(0);
  basisCache->setPhysicalCellNodes( _spectralConfusionMesh->physicalCellNodesForCell(0), cellIDs, true );

  double eps = .1; //
  FunctionPtr boundaryLayerFunction = Teuchos::rcp( new BoundaryLayerFunction(eps) );
  int numCells = basisCache->cellIDs().size();
  FieldContainer<double> integrals(numCells);
  double quadtol = 1e-2;
  double computedIntegral = boundaryLayerFunction->integrate(_spectralConfusionMesh,quadtol);
  double trueIntegral = (eps*(exp(1/eps) - exp(-1/eps)))*2.0;
  double diff = trueIntegral-computedIntegral;
  double relativeError = abs(diff)/abs(trueIntegral); // relative error

  double tol = 1e-2;
  if (relativeError > tol)
  {
    success = false;
    cout << "failing testAdaptiveIntegrate() with computed integral " << computedIntegral << " and true integral " << trueIntegral << endl;
  }
  return success;
}

  void values(FieldContainer<double> &values, BasisCachePtr basisCache){
    vector<int> cellIDs = basisCache->cellIDs();
    int numPoints = values.dimension(1);
    for (int i = 0;i<cellIDs.size();i++){
      double energyError = _energyErrorForCell[cellIDs[i]];
      for (int j = 0;j<numPoints;j++){
	values(i,j) = energyError;
      }
    }
  }

  void values(FieldContainer<double> &values, BasisCachePtr basisCache){
    vector<int> cellIDs = basisCache->cellIDs();
    int numPoints = values.dimension(1);
    for (int i = 0;i<cellIDs.size();i++){
      int partitionNumber = _mesh->partitionForCellID(cellIDs[i]);
      for (int j = 0;j<numPoints;j++){
	values(i,j) = partitionNumber;
      }
    }
  }

void MeshTransformationFunction::values(FieldContainer<double> &values, BasisCachePtr basisCache)
{
  CHECK_VALUES_RANK(values);
  vector<GlobalIndexType> cellIDs = basisCache->cellIDs();
  // identity map is the right thing most of the time
  // we'll do something different only where necessary
  int spaceDim = values.dimension(2);
  TEUCHOS_TEST_FOR_EXCEPTION(basisCache->cellTopology()->getDimension() != spaceDim, std::invalid_argument, "cellTopology dimension does not match the shape of the values container");
  if (_op == OP_VALUE)
  {
    values = basisCache->getPhysicalCubaturePoints(); // identity
  }
  else if (_op == OP_DX)
  {
    // identity map is 1 in all the x slots, 0 in all others
    int mod_value = 0; // the x slots are the mod spaceDim = 0 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
  else if (_op == OP_DY)
  {
    // identity map is 1 in all the y slots, 0 in all others
    int mod_value = 1; // the y slots are the mod spaceDim = 1 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
  else if (_op == OP_DZ)
  {
    // identity map is 1 in all the z slots, 0 in all others
    int mod_value = 2; // the z slots are the mod spaceDim = 2 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
//  if (_op == OP_DX) {
//    cout << "values before cellTransformation:\n" << values;
//  }
  for (int cellIndex=0; cellIndex < cellIDs.size(); cellIndex++)
  {
    GlobalIndexType cellID = cellIDs[cellIndex];
    if (_cellTransforms.find(cellID) == _cellTransforms.end()) continue;
    TFunctionPtr<double> cellTransformation = _cellTransforms[cellID];
    ((CellTransformationFunction*)cellTransformation.get())->setCellIndex(cellIndex);
    cellTransformation->values(values, basisCache);
  }
//  if (_op == OP_DX) {
//    cout << "values after cellTransformation:\n" << values;
//  }
}

void MeshPolyOrderFunction::values(FieldContainer<double> &values, BasisCachePtr basisCache) {
  vector<GlobalIndexType> cellIDs = basisCache->cellIDs();
  IndexType cellIndex = 0;
  int numPoints = values.dimension(1);
  for (vector<GlobalIndexType>::iterator cellIDIt = cellIDs.begin(); cellIDIt != cellIDs.end(); cellIDIt++, cellIndex++) {
    GlobalIndexType cellID = *cellIDIt;
    int polyOrder = _mesh->cellPolyOrder(cellID); // H1 order
    for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
      values(cellIndex, ptIndex) = polyOrder-1;
    }
  }
}

// computes riesz representation over a single element - map is from int (testID) to FieldContainer of values (sized cellIndex, numPoints)
void RieszRep::computeRepresentationValues(FieldContainer<double> &values, int testID, IntrepidExtendedTypes::EOperatorExtended op, BasisCachePtr basisCache){

  if (_repsNotComputed){
    cout << "Computing riesz rep dofs" << endl;
    computeRieszRep();
  }

  int spaceDim = _mesh->getTopology()->getSpaceDim();
  int numCells = values.dimension(0);
  int numPoints = values.dimension(1);
  vector<GlobalIndexType> cellIDs = basisCache->cellIDs();

  // all elems coming in should be of same type
  ElementPtr elem = _mesh->getElement(cellIDs[0]);
  ElementTypePtr elemTypePtr = elem->elementType();   
  DofOrderingPtr testOrderingPtr = elemTypePtr->testOrderPtr;
  CellTopoPtrLegacy cellTopoPtr = elemTypePtr->cellTopoPtr;
  int numTestDofsForVarID = testOrderingPtr->getBasisCardinality(testID, 0);
  BasisPtr testBasis = testOrderingPtr->getBasis(testID);
  
  bool testBasisIsVolumeBasis = (spaceDim == testBasis->domainTopology()->getDimension());  
  bool useCubPointsSideRefCell = testBasisIsVolumeBasis && basisCache->isSideCache();
  
  Teuchos::RCP< const FieldContainer<double> > transformedBasisValues = basisCache->getTransformedValues(testBasis,op,useCubPointsSideRefCell);
  
  int rank = values.rank() - 2; // if values are shaped as (C,P), scalar...
  if (rank > 1) {
    cout << "ranks greater than 1 not presently supported...\n";
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "ranks greater than 1 not presently supported...");
  }
  
//  Camellia::print("cellIDs",cellIDs);
  
  values.initialize(0.0);
  for (int cellIndex = 0;cellIndex<numCells;cellIndex++){
    int cellID = cellIDs[cellIndex];
    for (int j = 0;j<numTestDofsForVarID;j++) {
      int dofIndex = testOrderingPtr->getDofIndex(testID, j);
      for (int i = 0;i<numPoints;i++) {
        if (rank==0) {
          double basisValue = (*transformedBasisValues)(cellIndex,j,i);
          values(cellIndex,i) += basisValue*_rieszRepDofsGlobal[cellID](dofIndex);
        } else {
          for (int d = 0; d<spaceDim; d++) {
            double basisValue = (*transformedBasisValues)(cellIndex,j,i,d);
            values(cellIndex,i,d) += basisValue*_rieszRepDofsGlobal[cellID](dofIndex);
          }
        }
      }
    }
  }
//  TestSuite::serializeOutput("rep values", values);
}

 void values(FieldContainer<double> &values, BasisCachePtr basisCache)
 {
   vector<GlobalIndexType> cellIDs = basisCache->cellIDs();
   int numPoints = values.dimension(1);
   FieldContainer<double> points = basisCache->getPhysicalCubaturePoints();
   for (int i = 0; i<cellIDs.size(); i++)
   {
     for (int j = 0; j<numPoints; j++)
     {
       values(i,j) = cellIDs[i];
     }
   }
 }

void SideParityFunction::values(Intrepid::FieldContainer<double> &values, BasisCachePtr sideBasisCache)
{
  this->CHECK_VALUES_RANK(values);
  int numCells = values.dimension(0);
  int numPoints = values.dimension(1);
  int sideIndex = sideBasisCache->getSideIndex();
  if (sideIndex == -1)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "non-sideBasisCache passed into SideParityFunction");
  }
  if (sideBasisCache->getCellSideParities().size() > 0)
  {
    // then we'll use this, and won't require that mesh and cellIDs are set
    if (sideBasisCache->getCellSideParities().dimension(0) != numCells)
    {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "sideBasisCache->getCellSideParities() is non-empty, but the cell dimension doesn't match that of the values FieldContainer.");
    }

    for (int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++)
    {
      int parity = sideBasisCache->getCellSideParities()(cellOrdinal,sideIndex);
      for (int ptOrdinal=0; ptOrdinal<numPoints; ptOrdinal++)
      {
        values(cellOrdinal,ptOrdinal) = parity;
      }
    }
  }
  else
  {
    vector<GlobalIndexType> cellIDs = sideBasisCache->cellIDs();
    if (cellIDs.size() != numCells)
    {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "cellIDs.size() != numCells");
    }
    Teuchos::RCP<Mesh> mesh = sideBasisCache->mesh();
    if (! mesh.get())
    {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "mesh unset in BasisCache.");
    }
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      int parity = mesh->cellSideParitiesForCell(cellIDs[cellIndex])(0,sideIndex);
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        values(cellIndex,ptIndex) = parity;
      }
    }
  }
}

    virtual void localStiffnessMatrixAndRHS(FieldContainer<double> &localStiffness, FieldContainer<double> &rhsVector,
                                            IPPtr ip, BasisCachePtr ipBasisCache,
                                            RHSPtr rhs, BasisCachePtr basisCache) {
        double testMatrixAssemblyTime = 0, testMatrixInversionTime = 0, localStiffnessDeterminationFromTestsTime = 0;

#ifdef HAVE_MPI
        Epetra_MpiComm Comm(MPI_COMM_WORLD);
        //cout << "rank: " << rank << " of " << numProcs << endl;
#else
        Epetra_SerialComm Comm;
#endif

        Epetra_Time timer(Comm);

        // localStiffness should have dim. (numCells, numTrialFields, numTrialFields)
        MeshPtr mesh = basisCache->mesh();
        if (mesh.get() == NULL) {
            cout << "localStiffnessMatrix requires BasisCache to have mesh set.\n";
            TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "localStiffnessMatrix requires BasisCache to have mesh set.");
        }
        const vector<GlobalIndexType>* cellIDs = &basisCache->cellIDs();
        int numCells = cellIDs->size();
        if (numCells != localStiffness.dimension(0)) {
            cout << "localStiffnessMatrix requires basisCache->cellIDs() to have the same # of cells as the first dimension of localStiffness\n";
            TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "localStiffnessMatrix requires basisCache->cellIDs() to have the same # of cells as the first dimension of localStiffness");
        }

        ElementTypePtr elemType = mesh->getElementType((*cellIDs)[0]); // we assume all cells provided are of the same type
        DofOrderingPtr trialOrder = elemType->trialOrderPtr;
        DofOrderingPtr fieldOrder = mesh->getDofOrderingFactory().getFieldOrdering(trialOrder);
        DofOrderingPtr traceOrder = mesh->getDofOrderingFactory().getTraceOrdering(trialOrder);

        map<int, int> stiffnessIndexForTraceIndex;
        map<int, int> stiffnessIndexForFieldIndex;
        set<int> varIDs = trialOrder->getVarIDs();
        for (set<int>::iterator varIt = varIDs.begin(); varIt != varIDs.end(); varIt++) {
            int varID = *varIt;
            const vector<int>* sidesForVar = &trialOrder->getSidesForVarID(varID);
            bool isTrace = (sidesForVar->size() > 1);
            for (vector<int>::const_iterator sideIt = sidesForVar->begin(); sideIt != sidesForVar->end(); sideIt++) {
                int sideOrdinal = *sideIt;
                vector<int> dofIndices = trialOrder->getDofIndices(varID,sideOrdinal);
                if (isTrace) {
                    vector<int> traceDofIndices = traceOrder->getDofIndices(varID,sideOrdinal);
                    for (int i=0; i<traceDofIndices.size(); i++) {
                        stiffnessIndexForTraceIndex[traceDofIndices[i]] = dofIndices[i];
                    }
                } else {
                    vector<int> fieldDofIndices = fieldOrder->getDofIndices(varID);
                    for (int i=0; i<fieldDofIndices.size(); i++) {
                        stiffnessIndexForFieldIndex[fieldDofIndices[i]] = dofIndices[i];
                    }
                }
            }
        }

        int numTrialDofs = trialOrder->totalDofs();
        if ((numTrialDofs != localStiffness.dimension(1)) || (numTrialDofs != localStiffness.dimension(2))) {
            cout << "localStiffness should have dimensions (C,numTrialFields,numTrialFields).\n";
            TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "localStiffness should have dimensions (C,numTrialFields,numTrialFields).");
        }

        map<int,int> traceTestMap, fieldTestMap;
        int numEquations = _virtualTerms.getFieldTestVars().size();
        for (int eqn=0; eqn<numEquations; eqn++) {
            VarPtr testVar = _virtualTerms.getFieldTestVars()[eqn];
            VarPtr traceVar = _virtualTerms.getAssociatedTrace(testVar);
            VarPtr fieldVar = _virtualTerms.getAssociatedField(testVar);
            traceTestMap[traceVar->ID()] = testVar->ID();
            fieldTestMap[fieldVar->ID()] = testVar->ID();
        }

        int maxDegreeField = fieldOrder->maxBasisDegree();
        int testDegreeInterior = maxDegreeField + _virtualTerms.getTestEnrichment();
        int testDegreeTrace = testDegreeInterior + 2;

        cout << "ERROR in Virtual: getRelabeledDofOrdering() is commented out in DofOrderingFactory.  Need to rewrite for the new caching scheme.\n";
        DofOrderingPtr testOrderInterior; // = mesh->getDofOrderingFactory().getRelabeledDofOrdering(fieldOrder, fieldTestMap);
        testOrderInterior = mesh->getDofOrderingFactory().setBasisDegree(testOrderInterior, testDegreeInterior, false);
        DofOrderingPtr testOrderTrace = mesh->getDofOrderingFactory().setBasisDegree(testOrderInterior, testDegreeTrace, true); // this has a bunch of extra dofs (interior guys)

        map<int, int> remappedTraceIndices; // go from the index that includes the interior dofs to one that doesn't
        set<int> testIDs = testOrderTrace->getVarIDs();
        int testTraceIndex = 0;
        for (set<int>::iterator testIDIt = testIDs.begin(); testIDIt != testIDs.end(); testIDIt++) {
            int testID = *testIDIt;
            BasisPtr basis = testOrderTrace->getBasis(testID);
            vector<int> interiorDofs = basis->dofOrdinalsForInterior();
            for (int basisOrdinal=0; basisOrdinal<basis->getCardinality(); basisOrdinal++) {
                if (std::find(interiorDofs.begin(),interiorDofs.end(),basisOrdinal) == interiorDofs.end()) {
                    int dofIndex = testOrderTrace->getDofIndex(testID, basisOrdinal);
                    remappedTraceIndices[dofIndex] = testTraceIndex;
                    testTraceIndex++;
                }
            }
        }

//    DofOrderingPtr testOrderTrace = mesh->getDofOrderingFactory().getRelabeledDofOrdering(traceOrder, traceTestMap);
//    testOrderTrace = mesh->getDofOrderingFactory().setBasisDegree(testOrderTrace, testDegreeTrace);

//.........这里部分代码省略.........

void PenaltyMethodFilter::filter(FieldContainer<double> &localStiffnessMatrix, FieldContainer<double> &localRHSVector,
                                 BasisCachePtr basisCache, Teuchos::RCP<Mesh> mesh, Teuchos::RCP<BC> bc){ 
  
  // assumption: filter gets elements of all the same type  
  TEUCHOS_TEST_FOR_EXCEPTION(basisCache->cellIDs().size()==0,std::invalid_argument,"no cell IDs given to filter");
  
  ElementTypePtr elemTypePtr = mesh->getElement(basisCache->cellIDs()[0])->elementType();
  int numCells = localStiffnessMatrix.dimension(0);
  
  DofOrderingPtr trialOrderPtr = elemTypePtr->trialOrderPtr;
  
  unsigned numSides = CamelliaCellTools::getSideCount( *elemTypePtr->cellTopoPtr );
  // only allows for L2 inner products at the moment. 
  IntrepidExtendedTypes::EOperatorExtended trialOperator =  IntrepidExtendedTypes::OP_VALUE;
	
  // loop over sides first 
  for (unsigned int sideIndex = 0; sideIndex<numSides; sideIndex++){
    
    // GET INTEGRATION INFO - get cubature points and side normals to send to Constraints (Cell,Point, spaceDim)
    FieldContainer<double> sideCubPoints = basisCache->getPhysicalCubaturePointsForSide(sideIndex);
    FieldContainer<double> sideNormals = basisCache->getSideUnitNormals(sideIndex);        
    
    int numPts = sideCubPoints.dimension(1);
    
    // GET CONSTRAINT INFO
    vector<map<int, FieldContainer<double> > > constrCoeffsVector;
    vector<FieldContainer<double> > constraintValuesVector;
    vector<FieldContainer<bool> > imposeHereVector;
    _constraints->getConstraints(sideCubPoints,sideNormals,constrCoeffsVector,constraintValuesVector);
    
    //loop thru constraints
    int i = 0;
    for (vector<map<int,FieldContainer<double> > >::iterator constrIt = constrCoeffsVector.begin();
         constrIt !=constrCoeffsVector.end(); constrIt++) {
      map<int,FieldContainer<double> > constrCoeffs = *constrIt;
      FieldContainer<double> constrValues = constraintValuesVector[i];
      i++;
      
      double penaltyParameter = 1e7; // (single_precision)^(-1) - perhaps have this computed relative to terms in the matrix?
      
      for (map<int,FieldContainer<double> >::iterator constrTestIDIt = constrCoeffs.begin();
           constrTestIDIt !=constrCoeffs.end(); constrTestIDIt++) {
        pair<int,FieldContainer<double> > constrTestPair = *constrTestIDIt;
        int testTrialID = constrTestPair.first;
        
        // get basis to integrate for testing fxns
        BasisPtr testTrialBasis = trialOrderPtr->getBasis(testTrialID,sideIndex);
        FieldContainer<double> testTrialValuesTransformedWeighted = *(basisCache->getTransformedWeightedValues(testTrialBasis,trialOperator,
                                                                                                              sideIndex,false));
        // make copies b/c we can't fudge with return values from basisCache (const) - dimensions (Cell,Field - basis ordinal, Point)
        FieldContainer<double> testTrialValuesWeightedCopy = testTrialValuesTransformedWeighted;
        
        int numDofs2 = trialOrderPtr->getBasisCardinality(testTrialID,sideIndex); 
        for (int cellIndex=0; cellIndex<numCells; cellIndex++){
          for (int dofIndex=0; dofIndex<numDofs2; dofIndex++){
            for (int ptIndex=0; ptIndex<numPts; ptIndex++){
              testTrialValuesWeightedCopy(cellIndex, dofIndex, ptIndex) *= constrTestPair.second(cellIndex, ptIndex); // scale by constraint coeff
            }	   
          }
        }
        
        // loop thru pairs of trialIDs and constr coeffs
        for (map<int,FieldContainer<double> >::iterator constrIDIt = constrCoeffs.begin();
             constrIDIt !=constrCoeffs.end(); constrIDIt++) {
          pair<int,FieldContainer<double> > constrPair = *constrIDIt;
          int trialID = constrPair.first;
          
          // get basis to integrate
          BasisPtr trialBasis1 = trialOrderPtr->getBasis(trialID,sideIndex);
          // for trial: the value lives on the side, so we don't use the volume coords either:
          FieldContainer<double> trialValuesTransformed = *(basisCache->getTransformedValues(trialBasis1,trialOperator,sideIndex,false));
          // make copies b/c we can't fudge with return values from basisCache (const) - dimensions (Cell,Field - basis ordinal, Point)
          FieldContainer<double> trialValuesCopy = trialValuesTransformed;
          
          // transform trial values
          int numDofs1 = trialOrderPtr->getBasisCardinality(trialID,sideIndex); 
          for (int dofIndex=0; dofIndex<numDofs1; dofIndex++){
            for (int cellIndex=0; cellIndex<numCells; cellIndex++){
              for (int ptIndex=0; ptIndex<numPts; ptIndex++){
                trialValuesCopy(cellIndex, dofIndex, ptIndex) *= constrPair.second(cellIndex, ptIndex); // scale by constraint coeff
              }	   
            }
          }
          
          /////////////////////////////////////////////////////////////////////////////////////
          
          
          // integrate the transformed values, add them to the relevant trial/testTrialID dof combos
          FieldContainer<double> unweightedPenaltyMatrix(numCells,numDofs1,numDofs2);
          FunctionSpaceTools::integrate<double>(unweightedPenaltyMatrix,trialValuesCopy,testTrialValuesWeightedCopy,COMP_BLAS);
          
          for (int cellIndex=0; cellIndex<numCells; cellIndex++){
            for (int testDofIndex=0; testDofIndex<numDofs2; testDofIndex++){		
              int localTestDof = trialOrderPtr->getDofIndex(testTrialID, testDofIndex, sideIndex);
              for (int trialDofIndex=0; trialDofIndex<numDofs1; trialDofIndex++){		
                int localTrialDof = trialOrderPtr->getDofIndex(trialID, trialDofIndex, sideIndex);
                localStiffnessMatrix(cellIndex,localTrialDof,localTestDof) 
                                 += penaltyParameter*unweightedPenaltyMatrix(cellIndex,trialDofIndex,testDofIndex);
              }
            }
//.........这里部分代码省略.........