C++ SparseMatrix::RemoveRowsColumns方法代码示例

void MyFrame::ExportMassMatrix(bool fullMatrix)
{
  if ((!fullMatrix) && (!precomputationState.fixedVerticesAvailable))
  {
    this->errMsg( _T("Error"),  
      _T("No fixed vertices have been specified.") );
    return;
  }

  wxFileDialog *dlg = new wxFileDialog(this, _T("Export mass matrix"), uiState.currentWorkingDirectory, _T(""), _T("Mass Matrix Files(*.M)|*.M|All files(*.*)|*.*"), wxFD_SAVE /*| wxHIDE_READONLY*/, wxDefaultPosition);
  if ( dlg->ShowModal() == wxID_OK )
  {
    wxString massMatrixFilename( dlg->GetPath() );
    SaveCurrentWorkingDirectory(massMatrixFilename);
    if( !massMatrixFilename.empty() )
    {
      // create mass matrix

      SparseMatrix * massMatrix;
      GenerateMassMatrix::computeMassMatrix(precomputationState.simulationMesh, &massMatrix, true); 

      if (!fullMatrix)
      {
        // constrain the degrees of freedom
        int numConstrainedVertices = (int) (precomputationState.fixedVertices.size());
        int * constrainedDOFs = (int*) malloc (sizeof(int) * 3 * numConstrainedVertices);
        int i = 0;
        for(set<int> :: iterator iter = precomputationState.fixedVertices.begin(); 
          iter != precomputationState.fixedVertices.end(); iter++)
        {
          constrainedDOFs[3*i+0] = 3 * (*iter) + 1;
          constrainedDOFs[3*i+1] = 3 * (*iter) + 2;
          constrainedDOFs[3*i+2] = 3 * (*iter) + 3;
          i++;
        }

        int oneIndexed = 1;
        massMatrix->RemoveRowsColumns(
          3 * numConstrainedVertices, constrainedDOFs, oneIndexed);
        free(constrainedDOFs);
      }

      const char * filename = massMatrixFilename.mb_str();
      int code = massMatrix->Save((char*)filename);
      delete(massMatrix);
      if (code != 0)
      {
        this->errMsg( _T("Saving error"),  
          _T("Unable to save mass matrix to ") + massMatrixFilename );
        dlg->Destroy();
        return;
      }
    }
  }

  dlg->Destroy();
}

void MyFrame::ExportStiffnessMatrix(bool fullMatrix)
{
  if ((!fullMatrix) && (!precomputationState.fixedVerticesAvailable))
  {
    this->errMsg( _T("Error"),  
      _T("No fixed vertices have been specified.") );
    return;
  }

  wxFileDialog *dlg = new wxFileDialog(this, _T("Export stiffness matrix"), uiState.currentWorkingDirectory, _T(""), _T("Stiffness Matrix Files(*.K)|*.K|All files(*.*)|*.*"), wxFD_SAVE /*| wxHIDE_READONLY*/, wxDefaultPosition);
  if ( dlg->ShowModal() == wxID_OK )
  {
    wxString stiffnessMatrixFilename( dlg->GetPath() );
    SaveCurrentWorkingDirectory(stiffnessMatrixFilename);
    if( !stiffnessMatrixFilename.empty() )
    {
      // create stiffness matrix
      StVKElementABCD * precomputedIntegrals = StVKElementABCDLoader::load(precomputationState.simulationMesh);
      StVKInternalForces * internalForces = 
        new StVKInternalForces(precomputationState.simulationMesh, precomputedIntegrals);

      SparseMatrix * stiffnessMatrix;
      StVKStiffnessMatrix * stiffnessMatrixClass = new StVKStiffnessMatrix(internalForces);
      stiffnessMatrixClass->GetStiffnessMatrixTopology(&stiffnessMatrix);
      double * zero = (double*) calloc(3 * precomputationState.simulationMesh->getNumVertices(), sizeof(double));
      stiffnessMatrixClass->ComputeStiffnessMatrix(zero, stiffnessMatrix);

      free(zero);
      delete(precomputedIntegrals);
      delete(stiffnessMatrixClass);
      delete(internalForces);

      if (!fullMatrix)
      {
        // constrain the degrees of freedom
        int numConstrainedVertices = (int) (precomputationState.fixedVertices.size());
        int * constrainedDOFs = (int*) malloc (sizeof(int) * 3 * numConstrainedVertices);
        int i = 0;
        for(set<int> :: iterator iter = precomputationState.fixedVertices.begin(); 
          iter != precomputationState.fixedVertices.end(); iter++)
        {
          constrainedDOFs[3*i+0] = 3 * (*iter) + 1;
          constrainedDOFs[3*i+1] = 3 * (*iter) + 2;
          constrainedDOFs[3*i+2] = 3 * (*iter) + 3;
          i++;
        }

        int oneIndexed = 1;
        stiffnessMatrix->RemoveRowsColumns(
          3 * numConstrainedVertices, constrainedDOFs, oneIndexed);
        free(constrainedDOFs);
      }

      const char * filename = stiffnessMatrixFilename.mb_str();
      int code = stiffnessMatrix->Save((char*)filename);

      delete(stiffnessMatrix);

      if (code != 0)
      {
        this->errMsg( _T("Saving error"),  
          _T("Unable to save stiffness matrix to ") + stiffnessMatrixFilename );
        dlg->Destroy();
        return;
      }
    }
  }

  dlg->Destroy();
}

void * MyFrame::LinearModesWorker(
      int numDesiredModes,
      int * r, double ** frequencies_, double ** modes_ )
{
  *r = -1;

  // create mass matrix
  SparseMatrix * massMatrix;
  GenerateMassMatrix::computeMassMatrix(precomputationState.simulationMesh, &massMatrix, true);

  // create stiffness matrix
  StVKElementABCD * precomputedIntegrals = StVKElementABCDLoader::load(precomputationState.simulationMesh);
  StVKInternalForces * internalForces = 
    new StVKInternalForces(precomputationState.simulationMesh, precomputedIntegrals);

  SparseMatrix * stiffnessMatrix;
  StVKStiffnessMatrix * stiffnessMatrixClass = new StVKStiffnessMatrix(internalForces);
  stiffnessMatrixClass->GetStiffnessMatrixTopology(&stiffnessMatrix);
  double * zero = (double*) calloc(3 * precomputationState.simulationMesh->getNumVertices(), sizeof(double));
  stiffnessMatrixClass->ComputeStiffnessMatrix(zero, stiffnessMatrix);

  free(zero);
  delete(precomputedIntegrals);
  delete(stiffnessMatrixClass);
  delete(internalForces);

  // constrain the degrees of freedom
  int numConstrainedVertices = (int) (precomputationState.fixedVertices.size());
  int * constrainedDOFs = (int*) malloc (sizeof(int) * 3 * numConstrainedVertices);
  set<int> :: iterator iter;
  int i = 0;
  for(iter = precomputationState.fixedVertices.begin(); iter != precomputationState.fixedVertices.end(); iter++)
  {
    constrainedDOFs[3*i+0] = 3 * (*iter) + 1;
    constrainedDOFs[3*i+1] = 3 * (*iter) + 2;
    constrainedDOFs[3*i+2] = 3 * (*iter) + 3;
    i++;
  }

  int oneIndexed = 1;
  massMatrix->RemoveRowsColumns(
    3 * numConstrainedVertices, constrainedDOFs, oneIndexed);

  stiffnessMatrix->RemoveRowsColumns(
    3 * numConstrainedVertices, constrainedDOFs, oneIndexed);

  // call ARPACK

  double * frequenciesTemp = (double*) malloc (sizeof(double) * numDesiredModes);
  int numRetainedDOFs = stiffnessMatrix->Getn();
  double * modesTemp = (double*) malloc 
    (sizeof(double) * numDesiredModes * numRetainedDOFs);

  printf("Computing linear modes using ARPACK: ...\n");
  PerformanceCounter ARPACKCounter;
  double sigma = -1.0;

  int numLinearSolverThreads = wxThread::GetCPUCount();
  if (numLinearSolverThreads > 3)
    numLinearSolverThreads = 3; // diminished returns in solver beyond 3 threads

  //massMatrix->Save("MFactory");
  //stiffnessMatrix->Save("KFactory");

  ARPACKSolver generalizedEigenvalueProblem;
  int nconv = generalizedEigenvalueProblem.SolveGenEigShInv
    (stiffnessMatrix, massMatrix, 
     numDesiredModes, frequenciesTemp, 
     modesTemp, sigma, numLinearSolverThreads);

  ARPACKCounter.StopCounter();
  double ARPACKTime = ARPACKCounter.GetElapsedTime();
  printf("ARPACK time: %G s.\n", ARPACKTime); fflush(NULL);

  if (nconv < numDesiredModes)
  {
    free(modesTemp);
    free(frequenciesTemp);
    *r = -3;
    free(constrainedDOFs);
    delete(massMatrix);
    delete(stiffnessMatrix);
    return NULL;
  }

  int n3 = 3 * precomputationState.simulationMesh->getNumVertices();
  *frequencies_ = (double*) calloc (numDesiredModes, sizeof(double));
  *modes_ = (double*) calloc (numDesiredModes * n3, sizeof(double));

  for(int i=0; i<numDesiredModes; i++)
  {
    // insert zero rows into the computed modes
    int oneIndexed = 1;
    InsertRows(n3, &modesTemp[numRetainedDOFs*i], &((*modes_)[n3*i]), 
      3 * numConstrainedVertices, constrainedDOFs, oneIndexed);
  }

  for(int i=0; i<numDesiredModes; i++)
  {
    if (frequenciesTemp[i] <= 0)
//.........这里部分代码省略.........