Source

PetIGA / test / IGACreate.c

Lisandro Dalcin c723954 

Lisandro Dalcin c8fe37d 
































Lisandro Dalcin c723954 



Lisandro Dalcin 6cb9ffa 
Lisandro Dalcin c723954 
Lisandro Dalcin c8fe37d 
Lisandro Dalcin 1a5126a 
Lisandro Dalcin c8fe37d 
Lisandro Dalcin c723954 



Lisandro Dalcin 7827519 
Lisandro Dalcin c723954 







Lisandro Dalcin 7827519 

Lisandro Dalcin c723954 

Lisandro Dalcin c8fe37d 







Lisandro Dalcin 1a5126a 
Lisandro Dalcin c8fe37d 



Lisandro Dalcin 1a5126a 



Lisandro Dalcin c8fe37d 










Lisandro Dalcin 7827519 






Lisandro Dalcin 9e70147 










Lisandro Dalcin c723954 

Lisandro Dalcin 1a5126a 




Lisandro Dalcin c723954 

#include "petiga.h"

#undef  __FUNCT__
#define __FUNCT__ "Scalar"
PetscErrorCode Scalar(IGAPoint p,const PetscScalar U[],PetscInt n,PetscScalar *S,void *ctx)
{
  PetscInt i;
  for (i=0; i<n; i++) S[i] = (PetscScalar)1.0;
  return 0;
}


#undef  __FUNCT__
#define __FUNCT__ "System"
PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
{
  PetscInt dof = p->dof;
  PetscInt nen = p->nen;
  PetscReal *N = p->shape[0];
  PetscInt a,b,i,j;
  for (a=0; a<nen; a++) {
    PetscReal Na = N[a];
    for (b=0; b<nen; b++) {
      PetscReal Nb = N[b];
      for (i=0; i<dof; i++)
        for (j=0; j<dof; j++)
          if (i==j)
            K[a*dof*nen*dof+i*nen*dof+b*dof+j] = Na*Nb;
    }
    for (i=0; i<dof; i++)
      F[a*dof+i] = Na * 1;
  }
  return 0;
}

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char *argv[]) {

  PetscInt       dim,dof;
  IGA            iga;
  PetscScalar    s;
  PetscReal      xmin,xmax;
  Vec            b,x;
  Mat            A;
  KSP            ksp;
  SNES           snes;
  TS             ts;
  DM             dm;
  PetscErrorCode ierr;
  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);

  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);

  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
  if (dim < 1) {ierr = IGASetDim(iga,dim=3);CHKERRQ(ierr);}
  ierr = IGAGetDof(iga,&dof);CHKERRQ(ierr);
  if (dof < 1) {ierr = IGASetDof(iga,dof=1);CHKERRQ(ierr);}
  ierr = IGASetUp(iga);CHKERRQ(ierr);

  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
  ierr = IGAFormScalar(iga,x,1,&s,Scalar,0);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);

  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
  ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
  ierr = IGASetUserSystem(iga,System,0);CHKERRQ(ierr);
  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);;CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
  ierr = VecMin(x,0,&xmin);CHKERRQ(ierr);
  ierr = VecMax(x,0,&xmax);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);

  ierr = IGACreateSNES(iga,&snes);CHKERRQ(ierr);
  ierr = SNESDestroy(&snes);CHKERRQ(ierr);

  ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
  ierr = TSDestroy(&ts);CHKERRQ(ierr);

  ierr = IGACreateElemDM(iga,dof,&dm);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = IGACreateGeomDM(iga,dim,&dm);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = IGACreateNodeDM(iga,1,&dm);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);

  ierr = IGACreateWrapperDM(iga,&dm);CHKERRQ(ierr);
  ierr = DMIGAGetIGA(dm,&iga);CHKERRQ(ierr);
  ierr = DMIGASetIGA(dm,iga);CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(dm,&x);CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(dm,&b);CHKERRQ(ierr);
  ierr = DMCreateMatrix(dm,NULL,&A);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);

  ierr = IGADestroy(&iga);CHKERRQ(ierr);

  if ((xmax-xmin) > 1e-2) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Unexpected result: x_min=%G x_max=%G\n",
                       (PetscScalar)xmin,(PetscScalar)xmax);CHKERRQ(ierr);
  }

  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}