Commits

Lisandro Dalcin  committed 9116c8b

More API calls for IGAElement and IGAPoint

  • Participants
  • Parent commits 125cc16

Comments (0)

Files changed (4)

File include/petiga.h

 
   PetscReal *point;    /*   [nqp][dim]                */
   PetscReal *weight;   /*   [nqp]                     */
-  PetscReal *detJ;     /*   [nqp]                     */
+  PetscReal *detJac;   /*   [nqp]                     */
   PetscReal *jacobian; /*   [nqp][dim][dim]           */
   PetscReal *shape[4]; /*0: [nqp][nen]                */
                        /*1: [nqp][nen][dim]           */
 extern PetscErrorCode IGAElementEnd(IGAElement element);
 
 extern PetscErrorCode IGAElementBuildFix(IGAElement element);
-extern PetscErrorCode IGAElementBuildMap(IGAElement element);
+extern PetscErrorCode IGAElementBuildMapping(IGAElement element);
 extern PetscErrorCode IGAElementBuildQuadrature(IGAElement element);
 extern PetscErrorCode IGAElementBuildShapeFuns(IGAElement element);
 
 extern PetscErrorCode IGAElementGetIndex(IGAElement element,PetscInt *index);
-extern PetscErrorCode IGAElementGetInfo(IGAElement element,PetscInt *nqp,PetscInt *nen,PetscInt *dof,PetscInt *dim);
+extern PetscErrorCode IGAElementGetSizes(IGAElement element,PetscInt *nen,PetscInt *dof,PetscInt *nqp);
+extern PetscErrorCode IGAElementGetMapping(IGAElement element,PetscInt *nen,const PetscInt *mapping[]);
+extern PetscErrorCode IGAElementGetQuadrature(IGAElement element,PetscInt *nqp,PetscInt *dim,
+                                              const PetscReal *point[],const PetscReal *weigth[],
+                                              const PetscReal *detJac[]);
+extern PetscErrorCode IGAElementGetShapeFuns(IGAElement element,PetscInt *nqp,PetscInt *nen,PetscInt *dim,
+                                             const PetscReal *jacobian[],const PetscReal **shapefuns[]);
+
 extern PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point);
 
 extern PetscErrorCode IGAElementGetWorkVec(IGAElement element,PetscScalar *V[]);
 
   PetscReal *point;    /*   [dim] */
   PetscReal weight;    /*   []    */
-  PetscReal detJ;      /*   []    */
+  PetscReal detJac;    /*   []    */
   PetscReal *jacobian; /*   [dim][dim] */
   PetscReal *shape[4]; /*0: [nen]      */
                        /*1: [nen][dim] */
 extern PetscBool      IGAPointNext(IGAPoint point);
 
 extern PetscErrorCode IGAPointGetIndex(IGAPoint point,PetscInt *index);
-extern PetscErrorCode IGAPointGetInfo(IGAPoint point,PetscInt *nen,PetscInt *dof,PetscInt *dim);
+extern PetscErrorCode IGAPointGetSizes(IGAPoint point,PetscInt *nen,PetscInt *dof,PetscInt *dim);
+extern PetscErrorCode IGAPointGetQuadrature(IGAPoint pnt,PetscInt *dim,
+                                            const PetscReal *point[],PetscReal *weigth,PetscReal *detJac);
+extern PetscErrorCode IGAPointGetShapeFuns(IGAPoint pnt,PetscInt *nen,PetscInt *dim,
+                                           const PetscReal *jacobian[],const PetscReal **shapefuns[]);
 
 extern PetscErrorCode IGAPointGetWorkVec(IGAPoint point,PetscScalar *V[]);
 extern PetscErrorCode IGAPointGetWorkMat(IGAPoint point,PetscScalar *M[]);

File src/petigaelem.c

   ierr = PetscFree(element->geometry);CHKERRQ(ierr);
   ierr = PetscFree(element->point);CHKERRQ(ierr);
   ierr = PetscFree(element->weight);CHKERRQ(ierr);
-  ierr = PetscFree(element->detJ);CHKERRQ(ierr);
+  ierr = PetscFree(element->detJac);CHKERRQ(ierr);
   ierr = PetscFree(element->jacobian);CHKERRQ(ierr);
   ierr = PetscFree(element->shape[0]);CHKERRQ(ierr);
   ierr = PetscFree(element->shape[1]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*(dim+1),PetscReal,&element->geometry);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*dim,PetscReal,&element->point);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp,PetscReal,&element->weight);CHKERRQ(ierr);
-    ierr = PetscMalloc1(nqp,PetscReal,&element->detJ);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp,PetscReal,&element->detJac);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*dim*dim,PetscReal,&element->jacobian);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen,PetscReal,&element->shape[0]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen*dim,PetscReal,&element->shape[1]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen*dim*dim,PetscReal,&element->shape[2]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen*dim*dim*dim,PetscReal,&element->shape[3]);CHKERRQ(ierr);
 
-    ierr = PetscMemzero(element->detJ,    sizeof(PetscReal)*nqp);CHKERRQ(ierr);
+    ierr = PetscMemzero(element->detJac,  sizeof(PetscReal)*nqp);CHKERRQ(ierr);
     ierr = PetscMemzero(element->jacobian,sizeof(PetscReal)*nqp*dim*dim);CHKERRQ(ierr);
     ierr = PetscMemzero(element->shape[0],sizeof(PetscReal)*nqp*nen);CHKERRQ(ierr);
     ierr = PetscMemzero(element->shape[1],sizeof(PetscReal)*nqp*nen*dim);CHKERRQ(ierr);
     index = (index - coord) / size;
     ID[i] = coord + range[0];
   }
-  IGAElementBuildMap(element);
+  IGAElementBuildMapping(element);
   IGAElementBuildFix(element);
   return PETSC_TRUE;
 }
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementGetInfo"
-PetscErrorCode IGAElementGetInfo(IGAElement element,PetscInt *nqp, PetscInt *nen,PetscInt *dof,PetscInt *dim)
+#define __FUNCT__ "IGAElementGetSizes"
+PetscErrorCode IGAElementGetSizes(IGAElement element,PetscInt *nen,PetscInt *dof,PetscInt *nqp)
 {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
+  if (nen) PetscValidIntPointer(nen,2);
+  if (dof) PetscValidIntPointer(dof,3);
   if (nqp) PetscValidIntPointer(nqp,4);
-  if (nen) PetscValidIntPointer(nen,5);
-  if (dof) PetscValidIntPointer(dof,3);
-  if (dim) PetscValidIntPointer(dim,2);
-  if (nqp) *nqp = element->nqp;
   if (nen) *nen = element->nen;
   if (dof) *dof = element->dof;
-  if (dim) *dim = element->dim;
+  if (nqp) *nqp = element->nqp;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementGetMapping"
+PetscErrorCode IGAElementGetMapping(IGAElement element,PetscInt *nen,const PetscInt *mapping[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  if (nen)     PetscValidIntPointer(nen,3);
+  if (mapping) PetscValidPointer(mapping,3);
+  if (nen)     *nen     = element->nen;
+  if (mapping) *mapping = element->mapping;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementGetQuadrature"
+PetscErrorCode IGAElementGetQuadrature(IGAElement element,PetscInt *nqp,PetscInt *dim,
+                                       const PetscReal *point[],const PetscReal *weight[],
+                                       const PetscReal *detJac[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  if (nqp)    PetscValidIntPointer(nqp,2);
+  if (dim)    PetscValidIntPointer(dim,3);
+  if (point)  PetscValidPointer(point,4);
+  if (weight) PetscValidPointer(weight,5);
+  if (detJac) PetscValidPointer(detJac,6);
+  if (nqp)    *nqp    = element->nqp;
+  if (dim)    *dim    = element->dim;
+  if (point)  *point  = element->point;
+  if (weight) *weight = element->weight;
+  if (detJac) *detJac = element->detJac;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAElementGetShapeFuns"
+PetscErrorCode IGAElementGetShapeFuns(IGAElement element,PetscInt *nqp,PetscInt *nen,PetscInt *dim,
+                                      const PetscReal *jacobian[],const PetscReal **shapefuns[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  if (nqp)       PetscValidIntPointer(nqp,2);
+  if (nen)       PetscValidIntPointer(nen,3);
+  if (dim)       PetscValidIntPointer(dim,4);
+  if (jacobian)  PetscValidPointer(jacobian,5);
+  if (shapefuns) PetscValidPointer(shapefuns,6);
+  if (nqp)       *nqp       = element->nqp;
+  if (nen)       *nqp       = element->nqp;
+  if (dim)       *dim       = element->dim;
+  if (jacobian)  *jacobian  = element->jacobian;
+  if (shapefuns) *shapefuns = (const PetscReal **)element->shape;
   PetscFunctionReturn(0);
 }
 
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAElementBuildMap"
-PetscErrorCode IGAElementBuildMap(IGAElement element)
+#define __FUNCT__ "IGAElementBuildMapping"
+PetscErrorCode IGAElementBuildMapping(IGAElement element)
 {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
     IGABasis *BD  = element->parent->basis;
     PetscInt *ID  = element->ID;
     PetscReal *Cw = element->geometry;
-    PetscReal *dJ = element->detJ;
+    PetscReal *dJ = element->detJac;
     PetscReal *J  = element->jacobian;
     PetscReal **N = element->shape;
     switch (element->dim) {

File src/petigapc.c

     IGA          iga = 0;
     IGAElement   element;
     PetscInt     nen,dof;
-    PetscInt     n,*idx;
-    PetscScalar  *vals,*work,lwkopt;
+    PetscInt     n,*indices;
+    PetscScalar  *values,*work,lwkopt;
     PetscBLASInt m,*ipiv,info,lwork;
     PetscInt     start,end;
-    const PetscInt *lgmap;
+    const PetscInt *ltogmap;
+    const PetscInt *mapping;
     ISLocalToGlobalMapping map;
 
     ierr = PetscObjectQuery((PetscObject)A,"IGA",(PetscObject*)&iga);CHKERRQ(ierr);
     PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
 
     ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
-    ierr = IGAElementGetInfo(element,0,&nen,&dof,0);CHKERRQ(ierr);
+    ierr = IGAElementGetSizes(element,&nen,&dof,0);CHKERRQ(ierr);
 
     if (dof == 1) {
       ierr = MatGetLocalToGlobalMapping(A,&map,PETSC_NULL);CHKERRQ(ierr);
     } else {
       ierr = MatGetLocalToGlobalMappingBlock(A,&map,PETSC_NULL);CHKERRQ(ierr);
     }
-    ierr = ISLocalToGlobalMappingGetIndices(map,&lgmap);CHKERRQ(ierr);
+    ierr = ISLocalToGlobalMappingGetIndices(map,&ltogmap);CHKERRQ(ierr);
     ierr = MatGetOwnershipRange(A,&start,&end);CHKERRQ(ierr);
     start /= dof; end /= dof;
 
     n = nen*dof;
-    ierr = PetscMalloc1(  n,PetscInt,   &idx);CHKERRQ(ierr);
-    ierr = PetscMalloc1(n*n,PetscScalar,&vals);CHKERRQ(ierr);
+    ierr = PetscMalloc2(n,PetscInt,&indices,n*n,PetscScalar,&values);CHKERRQ(ierr);
     m = PetscBLASIntCast(n); lwork = -1; work = &lwkopt;
     ierr = PetscMalloc1(m,PetscBLASInt,&ipiv);CHKERRQ(ierr);
-    LAPACKgetri_(&m,vals,&m,ipiv,work,&lwork,&info);
+    LAPACKgetri_(&m,values,&m,ipiv,work,&lwork,&info);
     lwork = (info==0) ? (PetscBLASInt)work[0] : m*128;
     ierr = PetscMalloc1(lwork,PetscScalar,&work);CHKERRQ(ierr);
 
     ierr = MatZeroEntries(B);CHKERRQ(ierr);
     ierr = IGAElementBegin(element);CHKERRQ(ierr);
     while (IGAElementNext(element)) {
-      const PetscInt *indices = element->mapping;
-      m = n = ComputeOwnedGlobalIndices(lgmap,start,end,dof,nen,indices,idx);
+      ierr = IGAElementGetMapping(element,&nen,&mapping);CHKERRQ(ierr);
+      m = n = ComputeOwnedGlobalIndices(ltogmap,start,end,dof,nen,mapping,indices);
       /* get element matrix from global matrix */
-      ierr = MatGetValues(A,n,idx,n,idx,vals);;CHKERRQ(ierr);
+      ierr = MatGetValues(A,n,indices,n,indices,values);CHKERRQ(ierr);
       /* compute inverse of element matrix */
-      LAPACKgetrf_(&m,&m,vals,&m,ipiv,&info);
+      LAPACKgetrf_(&m,&m,values,&m,ipiv,&info);
       if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LAPACKgetrf_");
       if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero-pivot in LU factorization");
       ierr = PetscLogFlops((1/3.*n*n*n         +2/3.*n));CHKERRQ(ierr); /* multiplications */
       ierr = PetscLogFlops((1/3.*n*n*n-1/2.*n*n+1/6.*n));CHKERRQ(ierr); /* additions */
-      LAPACKgetri_(&m,vals,&m,ipiv,work,&lwork,&info);
+      LAPACKgetri_(&m,values,&m,ipiv,work,&lwork,&info);
       if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LAPACKgetri_");
       if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero-pivot in LU factorization");
       ierr = PetscLogFlops((2/3.*n*n*n+1/2.*n*n+5/6.*n));CHKERRQ(ierr); /* multiplications */
       ierr = PetscLogFlops((2/3.*n*n*n-3/2.*n*n+5/6.*n));CHKERRQ(ierr); /* additions */
       /* add values back into preconditioner matrix */
-      ierr = MatSetValues(B,n,idx,n,idx,vals,ADD_VALUES);CHKERRQ(ierr);
+      ierr = MatSetValues(B,n,indices,n,indices,values,ADD_VALUES);CHKERRQ(ierr);
       ierr = PetscLogFlops(n*n);CHKERRQ(ierr);
     }
     ierr = IGAElementEnd(element);CHKERRQ(ierr);
     ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
     ierr = MatAssemblyEnd  (B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 
-    ierr = ISLocalToGlobalMappingRestoreIndices(map,&lgmap);CHKERRQ(ierr);
-    ierr = PetscFree(idx);CHKERRQ(ierr);
-    ierr = PetscFree(vals);CHKERRQ(ierr);
+    ierr = ISLocalToGlobalMappingRestoreIndices(map,&ltogmap);CHKERRQ(ierr);
+    ierr = PetscFree2(indices,values);CHKERRQ(ierr);
     ierr = PetscFree(ipiv);CHKERRQ(ierr);
     ierr = PetscFree(work);CHKERRQ(ierr);
   }
   PetscFunctionBegin;
   ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
   if (!isascii) PetscFunctionReturn(0);
+  if (!ebe->mat) PetscFunctionReturn(0);
   ierr = PetscViewerASCIIPrintf(viewer,"element-by-element preconditioner matrix:\n");CHKERRQ(ierr);
   ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
   ierr = MatView(ebe->mat,viewer);CHKERRQ(ierr);

File src/petigapoint.c

   /* */
   point->point    = parent->point + index * dim;
   point->weight   = parent->weight[index];
-  point->detJ     = parent->detJ[index];
+  point->detJac   = parent->detJac[index];
   point->jacobian = parent->jacobian + index * dim*dim;
   point->shape[0] = parent->shape[0] + index * nen;
   point->shape[1] = parent->shape[1] + index * nen*dim;
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAPointGetInfo"
-PetscErrorCode IGAPointGetInfo(IGAPoint point,PetscInt *nen,PetscInt *dof,PetscInt *dim)
+#define __FUNCT__ "IGAPointGetSizes"
+PetscErrorCode IGAPointGetSizes(IGAPoint point,PetscInt *nen,PetscInt *dof,PetscInt *dim)
 {
   PetscFunctionBegin;
   PetscValidPointer(point,1);
   PetscFunctionReturn(0);
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "IGAPointGetQuadrature"
+PetscErrorCode IGAPointGetQuadrature(IGAPoint pnt,PetscInt *dim,
+                                     const PetscReal *point[],PetscReal *weight,PetscReal *detJac)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(pnt,1);
+  if (dim)    PetscValidIntPointer(dim,2);
+  if (point)  PetscValidPointer(point,3);
+  if (weight) PetscValidDoublePointer(weight,4);
+  if (detJac) PetscValidDoublePointer(detJac,5);
+  if (dim)    *dim    = pnt->dim;
+  if (point)  *point  = pnt->point;
+  if (weight) *weight = pnt->weight;
+  if (detJac) *detJac = pnt->detJac;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAPointGetShapeFuns"
+PetscErrorCode IGAPointGetShapeFuns(IGAPoint pnt,PetscInt *nen,PetscInt *dim,
+                                    const PetscReal *jacobian[],const PetscReal **shapefuns[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(pnt,1);
+  if (nen)       PetscValidIntPointer(nen,2);
+  if (dim)       PetscValidIntPointer(dim,3);
+  if (jacobian)  PetscValidPointer(jacobian,4);
+  if (shapefuns) PetscValidPointer(shapefuns,5);
+  if (nen)       *nen       = pnt->nen;
+  if (dim)       *dim       = pnt->dim;
+  if (jacobian)  *jacobian  = pnt->jacobian;
+  if (shapefuns) *shapefuns = (const PetscReal **)pnt->shape;
+  PetscFunctionReturn(0);
+}
+
 EXTERN_C_BEGIN
 extern void IGA_Interpolate(PetscInt nen,PetscInt dof,PetscInt dim,PetscInt der,
                             const PetscReal N[],const PetscScalar U[],PetscScalar u[]);
   {
     PetscInt nen = point->nen;
     PetscInt dof = point->dof;
-    PetscReal JW = point->detJ*point->weight;
+    PetscReal JW = point->detJac*point->weight;
     PetscInt i, n = nen*dof;
     for (i=0; i<n; i++) F[i] += f[i] * JW;
   }
   {
     PetscInt nen = point->nen;
     PetscInt dof = point->dof;
-    PetscReal JW = point->detJ*point->weight;
+    PetscReal JW = point->detJac*point->weight;
     PetscInt i, n = (nen*dof)*(nen*dof);
     for (i=0; i<n; i++) K[i] += k[i] * JW;
   }