Lisandro Dalcin avatar Lisandro Dalcin committed 57c32c9

Add IGASetFixTable()

Comments (0)

Files changed (6)

 PETSC_EXTERN PetscErrorCode IGAFormSetIEFunction(IGAForm form,IGAFormIEFunction IEFunction,void *ctx);
 PETSC_EXTERN PetscErrorCode IGAFormSetIEJacobian(IGAForm form,IGAFormIEJacobian IEJacobian,void *ctx);
 
+PETSC_EXTERN PetscErrorCode IGASetFixTable(IGA iga,Vec table);
 PETSC_EXTERN PetscErrorCode IGASetBoundaryValue(IGA iga,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
 PETSC_EXTERN PetscErrorCode IGASetBoundaryLoad (IGA iga,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
 PETSC_EXTERN PetscErrorCode IGASetBoundaryForm (IGA iga,PetscInt axis,PetscInt side,PetscBool flag);
   PetscBool   rational;
   PetscInt    geometry;
   PetscInt    property;
+  PetscBool   fixtable;
   PetscReal   *rationalW;
   PetscReal   *geometryX;
   PetscScalar *propertyA;
+  PetscScalar *fixtableU;
 
   PetscInt  proc_sizes[3];
   PetscInt  proc_ranks[3];
   iga->rational = PETSC_FALSE;
   iga->geometry = 0;
   iga->property = 0;
+  iga->fixtable = PETSC_FALSE;
   ierr = PetscFree(iga->rationalW);CHKERRQ(ierr);
   ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
   ierr = PetscFree(iga->propertyA);CHKERRQ(ierr);
+  ierr = PetscFree(iga->fixtableU);CHKERRQ(ierr);
   ierr = DMDestroy(&iga->geom_dm);CHKERRQ(ierr);
   /* node */
   ierr = AODestroy(&iga->ao);CHKERRQ(ierr);
   iga->rational = PETSC_FALSE;
   iga->geometry = 0;
   iga->property = 0;
+  iga->fixtable = PETSC_FALSE;
   ierr = PetscFree(iga->rationalW);CHKERRQ(ierr);
   ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
   ierr = PetscFree(iga->propertyA);CHKERRQ(ierr);
+  ierr = PetscFree(iga->fixtableU);CHKERRQ(ierr);
   ierr = DMDestroy(&iga->geom_dm);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 static void AddFixa(IGAElement element,IGAFormBC bc,PetscInt a)
 {
   if (bc->count) {
+    IGA iga = element->parent;
     PetscInt dof = element->dof;
     PetscInt count = element->nfix;
     PetscInt *index = element->ifix;
     PetscScalar *value = element->vfix;
     PetscInt j,k,n = bc->count;
     for (k=0; k<n; k++) {
-      PetscInt idx = a*dof + bc->field[k];
+      PetscInt c = bc->field[k];
+      PetscInt idx = a*dof + c;
       PetscScalar val = bc->value[k];
+      if (iga->fixtable) val = iga->fixtableU[c + element->mapping[a]*dof];
       for (j=0; j<count; j++)
         if (index[j] == idx) break;
       if (j == count) count++;
 /* --------------------------------------------------------------- */
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGASetFixTable"
+PetscErrorCode IGASetFixTable(IGA iga,Vec U)
+{
+  Vec               local;
+  PetscInt          nlocal;
+  const PetscScalar *vlocal;
+  PetscErrorCode    ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  if (iga->setupstage < 2) IGACheckSetUp(iga,1);
+
+  iga->fixtable = PETSC_FALSE;
+  ierr = PetscFree(iga->fixtableU);
+  if (!U) PetscFunctionReturn(0);
+
+  PetscValidHeaderSpecific(U,VEC_CLASSID,2);
+  ierr = IGAGetLocalVecArray(iga,U,&local,&vlocal);CHKERRQ(ierr);
+  ierr = VecGetSize(local,&nlocal);CHKERRQ(ierr);
+  iga->fixtable = PETSC_TRUE;
+  ierr = PetscMalloc1(nlocal,&iga->fixtableU);
+  ierr = PetscMemcpy(iga->fixtableU,vlocal,nlocal*sizeof(PetscScalar));CHKERRQ(ierr);
+  ierr = IGARestoreLocalVecArray(iga,U,&local,&vlocal);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGASetBoundaryValue"
 PetscErrorCode IGASetBoundaryValue(IGA iga,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value)
 {
+#include "petiga.h"
+
+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;
+}
+
+static PetscReal Solution(PetscInt dim,PetscReal x[])
+{
+  PetscInt i; PetscReal u = 0.0;
+  for (i=0; i<dim; i++) u += x[i]*x[i];
+  return u;
+}
+
+static PetscReal Forcing(PetscInt dim,PetscReal x[])
+{
+  PetscInt i; PetscReal f = 0.0;
+  for (i=0; i<dim; i++) f += -2.0;
+  return f;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "System1"
+PetscErrorCode System1(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt dim = p->dim;
+  PetscInt nen = p->nen;
+  PetscReal *N0 = (typeof(N0)) p->shape[0];
+
+  PetscReal x[3];
+  IGAPointFormGeomMap(p,x);
+  PetscReal g = Solution(dim,x);
+
+  PetscInt a,b;
+  for (a=0; a<nen; a++) {
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = N0[a]*N0[b];
+    F[a] = N0[a] * g;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "System2"
+PetscErrorCode System2(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt dim = p->dim;
+  PetscInt nen = p->nen;
+  PetscReal *N0        = (typeof(N0)) p->shape[0];
+  PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
+
+  PetscReal x[3];
+  IGAPointFormGeomMap(p,x);
+  PetscReal f = Forcing(dim,x);
+
+  PetscInt a,b;
+  for (a=0; a<nen; a++) {
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
+    F[a] = N0[a] * f;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+{
+  PetscScalar u;
+  IGAPointFormValue(p,U,&u);
+
+  PetscInt  dim = p->dim;
+  PetscReal x[3];
+  IGAPointFormGeomMap(p,x);
+  PetscReal g = Solution(dim,x);
+
+  PetscReal e = PetscAbsScalar(u - g);
+  S[0] = e*e;
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  PetscBool print_error = PETSC_FALSE;
+  PetscBool check_error = PETSC_FALSE;
+  PetscReal error_tol   = 1e-4;
+  PetscBool draw = PETSC_FALSE;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","FixTable Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-print_error","Prints the L2 error of the solution",__FILE__,print_error,&print_error,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-check_error","Checks the L2 error of the solution",__FILE__,error_tol,&error_tol,&check_error);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  if (iga->dim < 1) {ierr = IGASetDim(iga,2);CHKERRQ(ierr);}
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
+  ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
+
+  /* Solve L2 projection problem */
+  ierr = IGASetFormSystem(iga,System1,NULL);CHKERRQ(ierr);
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  {
+    KSP ksp;
+    ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
+    ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+    ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  }
+
+  /* Solve Poisson problem with Dirichlet BCs */
+  PetscInt dim,dir,side;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  for (dir=0; dir<dim; dir++)
+    for (side=0; side<2; side++)
+      {ierr = IGASetBoundaryValue(iga,dir,side,0,/*dummy*/0.0);CHKERRQ(ierr);}
+  ierr = IGASetFixTable(iga,x);CHKERRQ(ierr); /* Set vector to read BCs from */
+  ierr = IGASetFormSystem(iga,System2,NULL);CHKERRQ(ierr);
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  ierr = VecSet(x,0.0);CHKERRQ(ierr);
+  {
+    KSP ksp;
+    ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
+    ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+    ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  }
+  ierr = IGASetFixTable(iga,NULL);CHKERRQ(ierr); /* Clear vector to read BCs from */
+
+  PetscScalar error = 0;
+  ierr = IGAComputeScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
+  error = PetscSqrtReal(PetscRealPart(error));
+
+  if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);}
+  if (check_error) {if (PetscRealPart(error)>error_tol) SETERRQ2(PETSC_COMM_WORLD,1,"L2 error=%G > %G\n",error,error_tol);}
+  if (draw&&dim<3) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);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;
+}
 		 runex2a_1 runex2a_2 runex2a_3 runex2a_4 runex2a_8 \
 		 runex2a.rm IGAInputOutput.rm
 
+FixTable: FixTable.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
+runex4a_1:
+	-@${MPIEXEC} -n 1 ./FixTable ${OPTS} -iga_dim 3 -pc_type lu -check_error 1e-6 -iga_elements 1
+	-@${MPIEXEC} -n 1 ./FixTable ${OPTS} -iga_dim 1 -pc_type lu -check_error 1e-6
+runex4a_2:
+	-@${MPIEXEC} -n 1 ./FixTable ${OPTS} -iga_dim 3 -pc_type lu -check_error 1e-6 -iga_elements 1
+	-@${MPIEXEC} -n 1 ./FixTable ${OPTS} -iga_dim 2 -pc_type lu -check_error 1e-6
+runex4a_3:
+	-@${MPIEXEC} -n 1 ./FixTable ${OPTS} -iga_dim 3 -pc_type lu -check_error 1e-6 -iga_elements 1
+	-@${MPIEXEC} -n 1 ./FixTable ${OPTS} -iga_dim 3 -pc_type lu -check_error 1e-6 -iga_elements 4
+runex4b_1:
+	-@${MPIEXEC} -n 2 ./FixTable ${OPTS} -iga_dim 1 -ksp_rtol 1e-7 -check_error 1e-6
+	-@${MPIEXEC} -n 3 ./FixTable ${OPTS} -iga_dim 1 -ksp_rtol 1e-7 -check_error 1e-6
+runex4b_2:
+	-@${MPIEXEC} -n 4 ./FixTable ${OPTS} -iga_dim 2 -ksp_rtol 1e-7 -check_error 1e-6
+	-@${MPIEXEC} -n 9 ./FixTable ${OPTS} -iga_dim 2 -ksp_rtol 1e-7 -check_error 1e-6
+runex4b_3:
+	-@${MPIEXEC} -n  8 ./FixTable ${OPTS} -iga_dim 3 -iga_elements  8,8,8 -ksp_rtol 1e-7 -check_error 1e-6
+	-@${MPIEXEC} -n 12 ./FixTable ${OPTS} -iga_dim 3 -iga_elements 12,8,8 -ksp_rtol 1e-7 -check_error 1e-6
+FixTable = FixTable.PETSc \
+           runex4a_1 runex4a_2 runex4a_3 \
+           runex4b_1 runex4b_2 runex4b_3 \
+           FixTable.rm
 
 GeometryMap: GeometryMap.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 
 TESTEXAMPLES_C = $(IGACreate) \
 		 $(IGAInputOutput) \
+		 $(FixTable) \
 		 $(GeometryMap) \
 		 $(Test_SNES_2D) \
 		 $(Oscillator)
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.