Commits

Lisandro Dalcin committed 6c7f8ca

Add dimension-independet linear elasticity example

  • Participants
  • Parent commits e0da893

Comments (0)

Files changed (2)

demo/Elasticity.c

+/*
+  This code solves the 3D elasticity equations given the Lame
+  constants and subject to Dirichlet and Neumann boundary conditions.
+
+  keywords: steady, vector, linear
+ */
+#include "petiga.h"
+
+#if PETSC_VERSION_LT(3,5,0)
+#define KSPSetOperators(ksp,A,B) KSPSetOperators(ksp,A,B,SAME_NONZERO_PATTERN)
+#endif
+
+typedef struct {
+  PetscReal mu;
+  PetscReal lambda;
+} AppCtx;
+
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *KK,PetscScalar *FF,void *ctx)
+{
+  AppCtx    *user  = (AppCtx *)ctx;
+  PetscReal mu     = user->mu;
+  PetscReal lambda = user->lambda;
+
+  PetscInt a,b,nen = p->nen;
+  PetscInt i,j,dim = p->dim;
+
+  PetscReal   (*B)                = (typeof(B)) p->shape[0];
+  PetscReal   (*D)[dim]           = (typeof(D)) p->shape[1];
+  PetscScalar (*K)[dim][nen][dim] = (typeof(K)) KK;
+  PetscScalar (*F)[dim]           = (typeof(F)) FF;
+
+  for (a=0; a<nen; a++) {
+    for (b=0; b<nen; b++) {
+      PetscReal Kabii = DOT(dim,D[a],D[b]);
+      for (i=0; i<dim; i++)
+        K[a][i][b][i] += mu * Kabii;
+      for (i=0; i<dim; i++)
+        for (j=0; j<dim; j++)
+          K[a][i][b][j] += lambda * D[a][i]*D[b][j] + mu * D[a][j]*D[b][i];
+    }
+  }
+
+  for (a=0; a<nen; a++)
+    for (i=0; i<dim; i++)
+      F[a][i] = B[a] * 0.0;
+
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscInt       i,dim;
+  PetscErrorCode ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  // Define problem parameters
+  AppCtx user;
+  user.mu     = 1.0;
+  user.lambda = 1.0;
+  PetscInt  nld = 3;
+  PetscReal load[3] = {1.0, 0.0, 0.0};
+
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Elasticity Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsReal     ("-lambda","Lame's first parameter", __FILE__,user.lambda,&user.lambda,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal     ("-mu",    "Lame's second parameter",__FILE__,user.mu,    &user.mu,    NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsRealArray("-load",  "Boundary load",          __FILE__,load,&nld,               NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
+  ierr = IGAOptionsAlias("-d",  "2", "-iga_dim");
+  ierr = IGAOptionsAlias("-N", "32", "-iga_elements");
+  ierr = IGAOptionsAlias("-p", NULL, "-iga_degree");
+  ierr = IGAOptionsAlias("-k", NULL, "-iga_continuity");
+  ierr = IGAOptionsAlias("-q", NULL, "-iga_quadrature");
+
+  // Setup discretization
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,dim);CHKERRQ(ierr);
+  ierr = IGASetOrder(iga,1);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  // Set boundary conditions
+  for (i=0; i<dim; i++) {ierr = IGASetBoundaryValue(iga,0,0,i,0.0);CHKERRQ(ierr);}     // Dirichlet
+  for (i=0; i<dim; i++) {ierr = IGASetBoundaryLoad (iga,0,1,i,load[i]);CHKERRQ(ierr);} // Neumann
+
+  // Create linear system
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,System,&user);CHKERRQ(ierr);
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
+
+  // Solve linear system
+  KSP ksp;
+  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+
+  // Save geometry and solution vector
+  PetscBool save = PETSC_FALSE;
+  ierr = PetscOptionsGetBool(NULL,"-save",&save,NULL);CHKERRQ(ierr);
+  if (save) {
+    ierr = IGAWrite(iga,"Elasticity-geometry.dat");CHKERRQ(ierr);
+    ierr = IGAWriteVec(iga,x,"Elasticity-solution.dat");CHKERRQ(ierr);
+  }
+
+  // Draw solution vector
+  PetscBool draw = PETSC_FALSE;
+  PetscReal pause = 0.0;
+  ierr = PetscOptionsGetBool(NULL,"-draw",&draw,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsGetReal(NULL,"-pause",&pause,NULL);CHKERRQ(ierr);
+  if (draw && dim <= 2) {
+    ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
+    ierr = PetscSleep(pause);CHKERRQ(ierr);
+  }
+
+  // Cleanup
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = MatDestroy(&A);CHKERRQ(ierr);
+  ierr = VecDestroy(&x);CHKERRQ(ierr);
+  ierr = VecDestroy(&b);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}
 CahnHilliard3D \
 NavierStokesKorteweg2D \
 NavierStokesVMS \
+Elasticity \
 Elasticity3D \
 HyperElasticity \
 Richards \
 ElasticRod: ElasticRod.o ElasticRodFJ.o chkopts
 	${CLINKER} -o $@ $< ElasticRodFJ.o ${PETIGA_LIB}
 	${RM} -f $< ElasticRodFJ.o elasticrodfj.mod
+Elasticity: Elasticity.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
 Elasticity3D: Elasticity3D.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<