Lisandro Dalcin avatar Lisandro Dalcin committed fef0221

Add IGAComputeVector() and IGAComputeMatrix()

Comments (0)

Files changed (6)

 /* ---------------------------------------------------------------- */
 
 typedef PetscErrorCode (*IGAUserScalar)    (IGAPoint point,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx);
+typedef PetscErrorCode (*IGAUserVector)    (IGAPoint point,PetscScalar *F,void *ctx);
+typedef PetscErrorCode (*IGAUserMatrix)    (IGAPoint point,PetscScalar *K,void *ctx);
 typedef PetscErrorCode (*IGAUserSystem)    (IGAPoint point,PetscScalar *K,PetscScalar *F,void *ctx);
 typedef PetscErrorCode (*IGAUserFunction)  (IGAPoint point,const PetscScalar *U,PetscScalar *F,void *ctx);
 typedef PetscErrorCode (*IGAUserJacobian)  (IGAPoint point,const PetscScalar *U,PetscScalar *J,void *ctx);
 typedef struct _IGAUserOps *IGAUserOps;
 struct _IGAUserOps {
   /**/
+  IGAUserVector     Vector;
+  void              *VecCtx;
+  IGAUserMatrix     Matrix;
+  void              *MatCtx;
   IGAUserSystem     System;
   void              *SysCtx;
   /**/
 PETSC_EXTERN PetscErrorCode IGABoundarySetValue(IGABoundary boundary,PetscInt field,PetscScalar value);
 PETSC_EXTERN PetscErrorCode IGABoundarySetLoad (IGABoundary boundary,PetscInt field,PetscScalar value);
 
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserVector    (IGABoundary boundary,IGAUserVector     Vector,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGABoundarySetUserMatrix    (IGABoundary boundary,IGAUserMatrix     Matrix,    void *ctx);
 PETSC_EXTERN PetscErrorCode IGABoundarySetUserSystem    (IGABoundary boundary,IGAUserSystem     System,    void *ctx);
 PETSC_EXTERN PetscErrorCode IGABoundarySetUserFunction  (IGABoundary boundary,IGAUserFunction   Function,  void *ctx);
 PETSC_EXTERN PetscErrorCode IGABoundarySetUserJacobian  (IGABoundary boundary,IGAUserJacobian   Jacobian,  void *ctx);
 PETSC_EXTERN PetscErrorCode IGAGetLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);
 PETSC_EXTERN PetscErrorCode IGARestoreLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);
 
+PETSC_EXTERN PetscErrorCode IGASetUserVector    (IGA iga,IGAUserVector     Vector,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetUserMatrix    (IGA iga,IGAUserMatrix     Matrix,    void *ctx);
 PETSC_EXTERN PetscErrorCode IGASetUserSystem    (IGA iga,IGAUserSystem     System,    void *ctx);
 PETSC_EXTERN PetscErrorCode IGASetUserFunction  (IGA iga,IGAUserFunction   Function,  void *ctx);
 PETSC_EXTERN PetscErrorCode IGASetUserJacobian  (IGA iga,IGAUserJacobian   Jacobian,  void *ctx);
 /* ---------------------------------------------------------------- */
 
 PETSC_EXTERN PetscLogEvent IGA_FormScalar;
+PETSC_EXTERN PetscLogEvent IGA_FormVector;
+PETSC_EXTERN PetscLogEvent IGA_FormMatrix;
 PETSC_EXTERN PetscLogEvent IGA_FormSystem;
 PETSC_EXTERN PetscLogEvent IGA_FormFunction;
 PETSC_EXTERN PetscLogEvent IGA_FormJacobian;
 #define PCIGABBB "igabbb"
 
 PETSC_EXTERN PetscErrorCode IGACreateKSP(IGA iga,KSP *ksp);
+PETSC_EXTERN PetscErrorCode IGAComputeVector(IGA iga,Vec B);
+PETSC_EXTERN PetscErrorCode IGAComputeMatrix(IGA iga,Mat A);
 PETSC_EXTERN PetscErrorCode IGAComputeSystem(IGA iga,Mat A,Vec B);
 
 PETSC_EXTERN PetscErrorCode IGACreateSNES(IGA iga,SNES *snes);
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGASetUserVector"
+/*@
+   IGASetUserSystem - Set the user callback to form the vector
+   which represents the discretized L(w).
+
+   Logically collective on IGA
+
+   Input Parameters:
++  iga - the IGA context
+.  Vector - the function which evaluates L(w)
+-  ctx - user-defined context for evaluation routine (may be NULL)
+
+   Details of Vector:
+$  PetscErrorCode Vector(IGAPoint p,PetscScalar *F,void *ctx);
+
++  p - point at which to evaluate L(w)
+.  F - contribution to L(w)
+-  ctx - user-defined context for evaluation routine
+
+   Level: normal
+
+.keywords: IGA, setup linear system, vector assembly
+@*/
+PetscErrorCode IGASetUserVector(IGA iga,IGAUserVector Vector,void *VecCtx)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  if (Vector) iga->userops->Vector = Vector;
+  if (VecCtx) iga->userops->VecCtx = VecCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGASetUserMatrix"
+/*@
+   IGASetUserSystem - Set the user callback to form the matrix and vector
+   which represents the discretized a(w,u).
+
+   Logically collective on IGA
+
+   Input Parameters:
++  iga - the IGA context
+.  Matrix - the function which evaluates a(w,u)
+-  ctx - user-defined context for evaluation routine (may be NULL)
+
+   Details of System:
+$  PetscErrorCode System(IGAPoint p,PetscScalar *K,void *ctx);
+
++  p - point at which to evaluate a(w,u)
+.  K - contribution to a(w,u)
+-  ctx - user-defined context for evaluation routine
+
+   Level: normal
+
+.keywords: IGA, setup linear system, matrix assembly
+@*/
+PetscErrorCode IGASetUserMatrix(IGA iga,IGAUserMatrix Matrix,void *MatCtx)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  if (Matrix) iga->userops->Matrix = Matrix;
+  if (MatCtx) iga->userops->MatCtx = MatCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGASetUserSystem"
 /*@
    IGASetUserSystem - Set the user callback to form the matrix and vector

src/petigabound.c

     if (!(boundary)->userops) {_ierr = PetscNew(struct _IGAUserOps,&(boundary)->userops);CHKERRQ(_ierr);} \
   } while (0)
 
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserVector"
+PetscErrorCode IGABoundarySetUserVector(IGABoundary boundary,IGAUserVector Vector,void *VecCtx)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(boundary,1);
+  IGABoundaryEnsureUserOps(boundary);
+  if (Vector) boundary->userops->Vector = Vector;
+  if (VecCtx) boundary->userops->VecCtx = VecCtx;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetUserMatrix"
+PetscErrorCode IGABoundarySetUserMatrix(IGABoundary boundary,IGAUserMatrix Matrix,void *MatCtx)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(boundary,1);
+  IGABoundaryEnsureUserOps(boundary);
+  if (Matrix) boundary->userops->Matrix = Matrix;
+  if (MatCtx) boundary->userops->MatCtx = MatCtx;
+  PetscFunctionReturn(0);
+}
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGABoundarySetUserSystem"
 #include "petiga.h"
 
 PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserVector(IGAElement element,IGAUserVector *vec,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->Vector);
+  if (!ops) return PETSC_FALSE;
+  *vec = ops->Vector;
+  *ctx = ops->VecCtx;
+  return PETSC_TRUE;
+}
+
+PETSC_STATIC_INLINE
+PetscBool IGAElementNextUserMatrix(IGAElement element,IGAUserMatrix *mat,void **ctx)
+{
+  IGAUserOps ops;
+  while (IGAElementNextUserOps(element,&ops) && !ops->Matrix);
+  if (!ops) return PETSC_FALSE;
+  *mat = ops->Matrix;
+  *ctx = ops->MatCtx;
+  return PETSC_TRUE;
+}
+
+PETSC_STATIC_INLINE
 PetscBool IGAElementNextUserSystem(IGAElement element,IGAUserSystem *sys,void **ctx)
 {
   IGAUserOps ops;
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAComputeVector"
+PetscErrorCode IGAComputeVector(IGA iga,Vec vecB)
+{
+  IGAElement     element;
+  IGAPoint       point;
+  IGAUserVector  Vector;
+  void           *ctx;
+  PetscScalar    *B;
+  PetscScalar    *F;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(vecB,VEC_CLASSID,2);
+  IGACheckSetUp(iga,1);
+  IGACheckUserOp(iga,1,Vector);
+
+  ierr = VecZeroEntries(vecB);CHKERRQ(ierr);
+
+  ierr = PetscLogEventBegin(IGA_FormVector,iga,vecB,0,0);CHKERRQ(ierr);
+
+  /* Element loop */
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
+    ierr = IGAElementGetWorkVec(element,&B);CHKERRQ(ierr);
+    /* UserVector loop */
+    while (IGAElementNextUserVector(element,&Vector,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        ierr = IGAPointGetWorkVec(point,&F);CHKERRQ(ierr);
+        ierr = Vector(point,F,ctx);CHKERRQ(ierr);
+        ierr = IGAPointAddVec(point,F,B);CHKERRQ(ierr);
+      }
+      ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
+    }
+    ierr = IGAElementAssembleVec(element,B,vecB);CHKERRQ(ierr);
+  }
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
+
+  ierr = PetscLogEventEnd(IGA_FormVector,iga,vecB,0,0);CHKERRQ(ierr);
+
+  ierr = VecAssemblyBegin(vecB);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd  (vecB);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAComputeMatrix"
+PetscErrorCode IGAComputeMatrix(IGA iga,Mat matA)
+{
+  IGAElement     element;
+  IGAPoint       point;
+  IGAUserMatrix  Matrix;
+  void           *ctx;
+  PetscScalar    *A;
+  PetscScalar    *K;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(matA,MAT_CLASSID,2);
+  IGACheckSetUp(iga,1);
+  IGACheckUserOp(iga,1,Matrix);
+
+  ierr = MatZeroEntries(matA);CHKERRQ(ierr);
+
+  ierr = PetscLogEventBegin(IGA_FormMatrix,iga,matA,0,0);CHKERRQ(ierr);
+
+  /* Element loop */
+  ierr = IGABeginElement(iga,&element);CHKERRQ(ierr);
+  while (IGANextElement(iga,element)) {
+    ierr = IGAElementGetWorkMat(element,&A);CHKERRQ(ierr);
+    /* UserMatrix loop */
+    while (IGAElementNextUserMatrix(element,&Matrix,&ctx)) {
+      /* Quadrature loop */
+      ierr = IGAElementBeginPoint(element,&point);CHKERRQ(ierr);
+      while (IGAElementNextPoint(element,point)) {
+        ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
+        ierr = Matrix(point,K,ctx);CHKERRQ(ierr);
+        ierr = IGAPointAddMat(point,K,A);CHKERRQ(ierr);
+      }
+      ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
+    }
+    ierr = IGAElementAssembleMat(element,A,matA);CHKERRQ(ierr);
+  }
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
+
+  ierr = PetscLogEventEnd(IGA_FormMatrix,iga,matA,0,0);CHKERRQ(ierr);
+
+  ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd  (matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGAComputeSystem"
 /*@
    IGAComputeSystem - Form the matrix and vector which represents the
 PetscClassId  IGA_CLASSID = 0;
 
 PetscLogEvent IGA_FormScalar = 0;
+PetscLogEvent IGA_FormVector = 0;
+PetscLogEvent IGA_FormMatrix = 0;
 PetscLogEvent IGA_FormSystem = 0;
 PetscLogEvent IGA_FormFunction = 0;
 PetscLogEvent IGA_FormJacobian = 0;
   ierr = IGARegisterAll();CHKERRQ(ierr);
   /* Register Events */
   ierr = PetscLogEventRegister("IGAFormScalar",IGA_CLASSID,&IGA_FormScalar);CHKERRQ(ierr);
+  ierr = PetscLogEventRegister("IGAFormVector",IGA_CLASSID,&IGA_FormVector);CHKERRQ(ierr);
+  ierr = PetscLogEventRegister("IGAFormMatrix",IGA_CLASSID,&IGA_FormMatrix);CHKERRQ(ierr);
   ierr = PetscLogEventRegister("IGAFormSystem",IGA_CLASSID,&IGA_FormSystem);CHKERRQ(ierr);
   ierr = PetscLogEventRegister("IGAFormFunction",IGA_CLASSID,&IGA_FormFunction);CHKERRQ(ierr);
   ierr = PetscLogEventRegister("IGAFormJacobian",IGA_CLASSID,&IGA_FormJacobian);CHKERRQ(ierr);
   return 0;
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "Vector"
+PetscErrorCode Vector(IGAPoint p,PetscScalar *F,void *ctx)
+{
+  PetscInt dof = p->dof;
+  PetscInt nen = p->nen;
+  PetscReal *N = p->shape[0];
+  PetscInt a,i;
+  for (a=0; a<nen; a++) {
+    PetscReal Na = N[a];
+    for (i=0; i<dof; i++)
+      F[a*dof+i] = Na * 1;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Matrix"
+PetscErrorCode Matrix(IGAPoint p,PetscScalar *K,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;
+    }
+  }
+  return 0;
+}
 
 #undef  __FUNCT__
 #define __FUNCT__ "System"
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = KSPDestroy(&ksp);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 = IGASetUserVector(iga,Vector,0);CHKERRQ(ierr);
+  ierr = IGAComputeVector(iga,b);CHKERRQ(ierr);
+  ierr = IGASetUserMatrix(iga,Matrix,0);CHKERRQ(ierr);
+  ierr = IGAComputeMatrix(iga,A);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);
 
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.