Commits

Nathan Collier  committed 5e1094d

fixed bug in bcast, feature now functional for gradients also

  • Participants
  • Parent commits eec8d1b

Comments (0)

Files changed (2)

File src/petigapoint.c

    Input Parameters:
 +  iga - the IGA context
 .  U   - the vector
-.  p   - the parametric location in the IGA (size of iga->dim)
-.  u   - 
--  du  - 
+.  p   - the parametric location in the IGA [iga->dim]
+.  u   - pointer to the interpolated values [iga->dof]
+-  du  - pointer to the gradient of the interpolated values [iga->dof][iga->dim]
 
    Notes: This function is intended to be used in a monitor function
    if a few values of the solution are needed at specific
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidHeaderSpecific(U,VEC_CLASSID,2);
   PetscValidRealPointer(p,3);
-  PetscValidScalarPointer(u,4);
+  if(u)  PetscValidScalarPointer(u ,4);
+  if(du) PetscValidScalarPointer(du,5);
   PetscErrorCode     ierr;
   IGAElement         ele;
   IGAPoint           pnt;
   MPI_Comm           comm;
   PetscMPIInt        rank,lroot=0,root=0;
+  PetscMPIInt        lfound=0,found=0;
   Vec                localU;
   const PetscScalar *arrayU;
   PetscScalar       *Ul;
   ierr = IGAElementInit(ele,iga);CHKERRQ(ierr);
   pnt->point = p;
   pnt->dof   = iga->dof;
+  pnt->dim   = iga->dim;
  1. Nathan Collier author

    This was the problem. It only got defined when rank = root (inside the conditional) and thus I was broadcasting things of different sizes for the gradients.

 
   /* 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;
+    lroot  = rank;
+    lfound = 1;
     ierr = IGAPointInit(pnt,ele);CHKERRQ(ierr);
     ierr = IGAPointEval(iga,pnt);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVal(ele,&Ul);CHKERRQ(ierr);
     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(&lfound,&found,1,MPI_INT,MPI_SUM,comm);
+  if(found == 0){
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Point not found on any partition of the domain");
+  }
   ierr = MPI_Allreduce(&lroot,&root,1,MPI_INT,MPI_SUM,comm);
 
   /* broadcast result so all processors have same data */

File test/Test_Eval.c

 #include "petiga.h"
+#include "time.h"
 
 #define SQ(A) ((A)*(A))
 
   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;
+
+  PetscReal  pnt[2];
+  PetscScalar sol=0,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[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);