C++ NonlinearImplicitSystem::get_mesh方法代码示例

static PetscErrorCode  DMCreateDomainDecomposition_libMesh(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
{
  PetscFunctionBegin;
  PetscErrorCode ierr;
  DM_libMesh     *dlm = (DM_libMesh *)(dm->data);
  NonlinearImplicitSystem *sys = dlm->sys;
  IS emb;
  if(dlm->decomposition_type != DMLIBMESH_DOMAIN_DECOMPOSITION) PetscFunctionReturn(0);
  *len = dlm->decomposition->size();
  if(namelist)      {ierr = PetscMalloc(*len*sizeof(char*), namelist);  CHKERRQ(ierr);}
  if(innerislist)   {ierr = PetscMalloc(*len*sizeof(IS),    innerislist);    CHKERRQ(ierr);}
  if(outerislist)   *outerislist = PETSC_NULL; /* FIX: allow mesh-based overlap. */
  if(dmlist)        {ierr = PetscMalloc(*len*sizeof(DM),    dmlist);    CHKERRQ(ierr);}
  for(unsigned int d = 0; d < dlm->decomposition->size(); ++d) {
    std::set<numeric_index_type>               dindices;
    std::string                          dname;
    std::map<std::string, unsigned int>  dblockids;
    std::map<unsigned int,std::string>   dblocknames;
    unsigned int                         dbcount = 0;
    for(std::set<unsigned int>::const_iterator bit = (*dlm->decomposition)[d].begin(); bit != (*dlm->decomposition)[d].end(); ++bit){
      unsigned int b = *bit;
      std::string bname = (*dlm->blocknames)[b];
      dblockids.insert(std::pair<std::string, unsigned int>(bname,b));
      dblocknames.insert(std::pair<unsigned int,std::string>(b,bname));
      if(!dbcount) dname = bname;
      else   dname += "_" + bname;
      ++dbcount;
      if(!innerislist) continue;
      MeshBase::const_element_iterator       el     = sys->get_mesh().active_local_subdomain_elements_begin(b);
      const MeshBase::const_element_iterator end_el = sys->get_mesh().active_local_subdomain_elements_end(b);
      for ( ; el != end_el; ++el) {
        const Elem* elem = *el;
        std::vector<numeric_index_type> evindices;
        /* Iterate only over this DM's variables. */
        for(std::map<std::string, unsigned int>::const_iterator vit = dlm->varids->begin(); vit != dlm->varids->end(); ++vit) {
          unsigned int v = vit->second;
          // Get the degree of freedom indices for the given variable off the current element.
          sys->get_dof_map().dof_indices(elem, evindices, v);
          for(unsigned int i = 0; i < evindices.size(); ++i) {
            numeric_index_type dof = evindices[i];
            if(dof >= sys->get_dof_map().first_dof() && dof < sys->get_dof_map().end_dof()) /* might want to use variable_first/last_local_dof instead */
              dindices.insert(dof);
          }
        }
      }
    }
    if(namelist) {
      ierr = PetscStrallocpy(dname.c_str(),(*namelist)+d);            CHKERRQ(ierr);
    }
    if(innerislist) {
      PetscInt *darray;
      IS dis;
      ierr = PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray); CHKERRQ(ierr);
      numeric_index_type i = 0;
      for(std::set<numeric_index_type>::const_iterator it = dindices.begin(); it != dindices.end(); ++it) {
        darray[i] = *it;
        ++i;
      }
      ierr = ISCreateGeneral(((PetscObject)dm)->comm, dindices.size(),darray, PETSC_OWN_POINTER, &dis); CHKERRQ(ierr);
      if(dlm->embedding) {
        /* Create a relative embedding into the parent's index space. */
#if PETSC_RELEASE_LESS_THAN(3,3,1)
        ierr = ISMapFactorRight(dis,dlm->embedding, PETSC_TRUE, &emb); CHKERRQ(ierr);
#else
        ierr = ISEmbed(dis,dlm->embedding, PETSC_TRUE, &emb); CHKERRQ(ierr);
#endif
        PetscInt elen, dlen;
        ierr = ISGetLocalSize(emb, &elen); CHKERRQ(ierr);
        ierr = ISGetLocalSize(dis, &dlen);  CHKERRQ(ierr);
        if(elen != dlen) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed field %D", d);
        ierr = ISDestroy(&dis); CHKERRQ(ierr);
        dis = emb;
      }
      else {
        emb = dis;
      }
      if(innerislist) {
        ierr = PetscObjectReference((PetscObject)dis); CHKERRQ(ierr);
        (*innerislist)[d] = dis;
      }
      ierr = ISDestroy(&dis); CHKERRQ(ierr);
    }
    if(dmlist) {
      DM ddm;
      ierr = DMCreate(((PetscObject)dm)->comm, &ddm); CHKERRQ(ierr);
      ierr = DMSetType(ddm, DMLIBMESH);               CHKERRQ(ierr);
      DM_libMesh *ddlm = (DM_libMesh*)(ddm->data);
      ddlm->sys = dlm->sys;
      /* copy over the varids and varnames */
      *ddlm->varids    = *dlm->varids;
      *ddlm->varnames  = *dlm->varnames;
      /* set the blocks from the d-th part of the decomposition. */
      *ddlm->blockids    = dblockids;
      *ddlm->blocknames  = dblocknames;
      ierr = PetscObjectReference((PetscObject)emb); CHKERRQ(ierr);
      ddlm->embedding = emb;
      ddlm->embedding_type = DMLIBMESH_DOMAIN_EMBEDDING;

      ierr = DMLibMeshSetUpName_Private(ddm); CHKERRQ(ierr);
      ierr = DMSetFromOptions(ddm);           CHKERRQ(ierr);
//.........这里部分代码省略.........

static PetscErrorCode  DMCreateFieldDecomposition_libMesh(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
{
  PetscFunctionBegin;
  PetscErrorCode ierr;
  DM_libMesh     *dlm = (DM_libMesh *)(dm->data);
  NonlinearImplicitSystem *sys = dlm->sys;
  IS emb;
  if(dlm->decomposition_type != DMLIBMESH_FIELD_DECOMPOSITION) PetscFunctionReturn(0);

  *len = dlm->decomposition->size();
  if(namelist) {ierr = PetscMalloc(*len*sizeof(char*), namelist);  CHKERRQ(ierr);}
  if(islist)   {ierr = PetscMalloc(*len*sizeof(IS),    islist);    CHKERRQ(ierr);}
  if(dmlist)   {ierr = PetscMalloc(*len*sizeof(DM),    dmlist);    CHKERRQ(ierr);}
  DofMap& dofmap = dlm->sys->get_dof_map();
  for(unsigned int d = 0; d < dlm->decomposition->size(); ++d) {
    std::set<unsigned int>               dindices;
    std::string                          dname;
    std::map<std::string, unsigned int>  dvarids;
    std::map<unsigned int, std::string>  dvarnames;
    unsigned int                         dvcount = 0;
    for(std::set<unsigned int>::const_iterator dvit = (*dlm->decomposition)[d].begin(); dvit != (*dlm->decomposition)[d].end(); ++dvit){
      unsigned int v = *dvit;
      std::string vname = (*dlm->varnames)[v];
      dvarids.insert(std::pair<std::string, unsigned int>(vname,v));
      dvarnames.insert(std::pair<unsigned int,std::string>(v,vname));
      if(!dvcount) dname = vname;
      else   dname += "_" + vname;
      ++dvcount;
      if(!islist) continue;
      /* Iterate only over this DM's blocks. */
      for(std::map<std::string, unsigned int>::const_iterator bit = dlm->blockids->begin(); bit != dlm->blockids->end(); ++bit) {
	unsigned int b = bit->second;
	MeshBase::const_element_iterator el     = sys->get_mesh().active_local_subdomain_elements_begin(b);
	MeshBase::const_element_iterator end_el = sys->get_mesh().active_local_subdomain_elements_end(b);
	for ( ; el != end_el; ++el) {
	  const Elem* elem = *el;
	  //unsigned int e_subdomain = elem->subdomain_id();
	  std::vector<unsigned int> evindices;
	  // Get the degree of freedom indices for the given variable off the current element.  
	  dofmap.dof_indices(elem, evindices, v);
	  for(unsigned int i = 0; i < evindices.size(); ++i) {
	    unsigned int dof = evindices[i];
	    if(dof >= dofmap.first_dof() && dof < dofmap.end_dof()) /* might want to use variable_first/last_local_dof instead */
	      dindices.insert(dof);
	  }
	}
      }
    }
    if(namelist) {
      ierr = PetscStrallocpy(dname.c_str(),(*namelist)+d);            CHKERRQ(ierr);
    }
    if(islist) {
      IS dis;
      PetscInt *darray;
      ierr = PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray); CHKERRQ(ierr);
      unsigned int i = 0;
      for(std::set<unsigned int>::const_iterator it = dindices.begin(); it != dindices.end(); ++it) {
	darray[i] = *it;
	++i;
      }
      ierr = ISCreateGeneral(((PetscObject)dm)->comm, dindices.size(),darray, PETSC_OWN_POINTER, &dis); CHKERRQ(ierr);
      if(dlm->embedding) {
	/* Create a relative embedding into the parent's index space. */
	ierr = ISMapFactorRight(dis,dlm->embedding, PETSC_TRUE, &emb); CHKERRQ(ierr);
	PetscInt elen, dlen;
	ierr = ISGetLocalSize(emb, &elen); CHKERRQ(ierr);
	ierr = ISGetLocalSize(dis, &dlen); CHKERRQ(ierr);
	if(elen != dlen) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed subdomain %D", d);
	ierr = ISDestroy(&dis); CHKERRQ(ierr);
	dis = emb;
      }
      else {
	emb = dis;
      }
      (*islist)[d] = dis;
    }
    if(dmlist) {
      DM ddm;
      ierr = DMCreate(((PetscObject)dm)->comm, &ddm); CHKERRQ(ierr);
      ierr = DMSetType(ddm, DMLIBMESH);               CHKERRQ(ierr);
      DM_libMesh *ddlm = (DM_libMesh*)(ddm->data);
      ddlm->sys = dlm->sys;
      /* copy over the block ids and names */
      *ddlm->blockids = *dlm->blockids; 
      *ddlm->blocknames = *dlm->blocknames;
      /* set the vars from the d-th part of the decomposition. */
      *ddlm->varids     = dvarids;        
      *ddlm->varnames   = dvarnames;        
      ierr = PetscObjectReference((PetscObject)emb); CHKERRQ(ierr);
      ddlm->embedding = emb;
      ddlm->embedding_type = DMLIBMESH_FIELD_EMBEDDING;

      ierr = DMLibMeshSetUpName_Private(ddm); CHKERRQ(ierr);
      ierr = DMSetFromOptions(ddm);           CHKERRQ(ierr);
      (*dmlist)[d] = ddm;
    }
  }
  PetscFunctionReturn(0);
}