Nathan Collier avatar Nathan Collier committed 4689564

added interpolate functionality

Comments (0)

Files changed (4)

 PETSC_EXTERN PetscErrorCode IGASetUserIEFunction(IGA iga,IGAUserIEFunction IEFunction,void *ctx);
 PETSC_EXTERN PetscErrorCode IGASetUserIEJacobian(IGA iga,IGAUserIEJacobian IEJacobian,void *ctx);
 
+PETSC_EXTERN PetscBool IGALocateElement(IGA iga,PetscReal *pnt,IGAElement element);
+PETSC_EXTERN PetscErrorCode IGAPointEval(IGA iga,IGAPoint point);
+PETSC_EXTERN PetscErrorCode IGAInterpolate(IGA iga,Vec U,PetscReal p[],PetscScalar u[],PetscScalar du[]);
+
 #undef  DMIGA
 #define DMIGA "iga"
 PETSC_EXTERN PetscErrorCode IGACreateWrapperDM(IGA iga,DM *dm);
 
 /* ---------------------------------------------------------------- */
 
-PETSC_EXTERN PetscBool IGALocateElement(IGA iga,PetscReal *pnt,IGAElement element);
-PETSC_EXTERN PetscErrorCode IGAPointEval(IGA iga,IGAPoint point);
-
-/* ---------------------------------------------------------------- */
-
 #if defined(PETSC_USE_DEBUG)
 #define IGACheckSetUp(iga,arg) do {                                      \
     if (PetscUnlikely(!(iga)->setup))                                    \

src/petigapoint.c

 
 .keywords: IGA, locate, point, element
 @*/
-PetscBool IGALocateElement(IGA iga,PetscScalar *pnt,IGAElement element)
+PetscBool IGALocateElement(IGA iga,PetscReal *pnt,IGAElement element)
 {
   PetscErrorCode ierr;
   PetscInt i,j,e,m,dim=iga->dim,*ID = element->ID;
 +  iga - the IGA context
 -  point - the IGAPoint context
 
-   Notes:
-   The point assumes you have already called IGALocateElement and it
-   returned true. This function is intended to be used when
-   integrating on a manifold embedded in the IGA domain.
+   Notes: The point assumes that the point is contained in an element
+   on this processors partition (call IGALocateElement)
 
    Level: devel
 
   for(i=0;i<element->dim;i++){ ierr = PetscFree(basis[i]);CHKERRQ(ierr); }
   PetscFunctionReturn(0);
 }
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAInterpolate"
+/*@
+   IGAInterpolate - Interpolates a vector given the IGA at the
+   requested point in the domain.
+
+   Input Parameters:
++  iga - the IGA context
+.  U   - the vector
+.  p   - the parametric location in the IGA (size of iga->dim)
+.  u   - 
+-  du  - 
+
+   Notes: This function is intended to be used in a monitor function
+   if a few values of the solution are needed at specific
+   locations. Using it for interpolating a large number of points will
+   be slow as it needs to create and destroy several objects
+   internally.
+
+   Level: normal
+
+.keywords: interpolate
+@*/
+PetscErrorCode IGAInterpolate(IGA iga,Vec U,PetscReal p[],PetscScalar u[],PetscScalar du[])
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(U,VEC_CLASSID,2);
+  PetscValidRealPointer(p,3);
+  PetscValidScalarPointer(u,4);
+  PetscErrorCode     ierr;
+  IGAElement         ele;
+  IGAPoint           pnt;
+  MPI_Comm           comm;
+  PetscMPIInt        rank,lroot=0,root=0;
+  Vec                localU;
+  const PetscScalar *arrayU;
+  PetscScalar       *Ul;
+
+  /* get comm and rank */
+  ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
+  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
+
+  /* create the IGAPoint and IGAElement */
+  ierr = IGAPointCreate(&pnt);CHKERRQ(ierr);
+  ierr = IGAElementCreate(&ele);CHKERRQ(ierr);
+  ierr = IGAElementInit(ele,iga);CHKERRQ(ierr);
+  pnt->point = p;
+  pnt->dof   = iga->dof;
+
+  /* if point is on this processor, this rank is root */
+  ierr = IGAGetLocalVecArray(iga,U,&localU,&arrayU);CHKERRQ(ierr);
+  if(IGALocateElement(iga,pnt->point,ele)){
+    lroot = rank;
+    ierr = IGAPointInit(pnt,ele);CHKERRQ(ierr);
+    ierr = IGAPointEval(iga,pnt);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(ele,&Ul);CHKERRQ(ierr);
+    ierr = IGAElementGetValues(ele,arrayU,Ul);CHKERRQ(ierr);
+    if(u)  IGA_GetValue(pnt->nen,pnt->dof,pnt->shape[0],Ul,u);
+    if(du) IGA_GetGrad(pnt->nen,pnt->dof,pnt->dim,pnt->shape[1],Ul,du);
+  }
+  ierr = IGARestoreLocalVecArray(iga,U,&localU,&arrayU);CHKERRQ(ierr);
+  ierr = MPI_Allreduce(&lroot,&root,1,MPI_INT,MPI_SUM,comm);
+
+  /* broadcast result so all processors have same data */
+  if(u)  { ierr = MPI_Bcast( u,pnt->dof,         MPIU_SCALAR,root,comm);CHKERRQ(ierr); }
+  if(du) { ierr = MPI_Bcast(du,pnt->dof*pnt->dim,MPIU_SCALAR,root,comm);CHKERRQ(ierr); }
+  ierr = IGAElementDestroy(&ele);CHKERRQ(ierr);
+  ierr = IGAPointDestroy(&pnt);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+#include "petiga.h"
+
+#define SQ(A) ((A)*(A))
+
+PetscScalar CrossTerm(PetscReal x, PetscReal y)
+{
+  return x*y;
+}
+
+typedef struct {
+  PetscScalar (*Function)(PetscReal x, PetscReal y);
+} AppCtx;
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal x = p->point[0];
+  PetscReal y = p->point[1];
+  PetscScalar f = user->Function(x,y);
+
+  PetscReal *N = p->shape[0];
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na = N[a];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb = N[b];
+      K[a*nen+b] = Na * Nb;
+    }
+    F[a] = Na * f;
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  AppCtx user;
+  user.Function = CrossTerm;
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,2);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);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 = IGASetUserSystem(iga,System,&user);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);
+  
+  PetscScalar pnt[2],sol,gsol[2],err,gerr;
+  srand(time(NULL));
+  PetscInt i;
+  for(i=0;i<10;i++){
+    pnt[0] = rand();
+    pnt[1] = rand();
+    pnt[0] /= RAND_MAX; 
+    pnt[1] /= RAND_MAX;
+    PetscPrintf(PETSC_COMM_WORLD,"Evaluating solution at x = (%f,%f)\n",pnt[0],pnt[1]);
+    ierr = IGAInterpolate(iga,x,pnt,&sol,gsol);
+    err = fabs(sol-user.Function(pnt[0],pnt[1]));
+    gerr = sqrt(SQ(gsol[0]-pnt[1])+SQ(gsol[1]-pnt[0]));
+    PetscPrintf(PETSC_COMM_WORLD,"\t     Solution  =  %e                Error = %e\n",sol,err);
+    PetscPrintf(PETSC_COMM_WORLD,"\tgrad(Solution) = [%e %e]  Error = %e\n",gsol[0],gsol[1],gerr);
+  }
+
+  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;
+}
 
 OPTS=-nox -malloc_debug -malloc_dump
 
+Test_Eval: Test_Eval.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
 
 IGACreate: IGACreate.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
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.