Commits

Lisandro Dalcin committed d67d90c

Update handling of geometry and property

Comments (0)

Files changed (4)

   element->atboundary  = PETSC_FALSE;
   element->boundary_id = -1;
 
+  element->rational = iga->rational ? PETSC_TRUE : PETSC_FALSE;
+  element->geometry = iga->geometry ? PETSC_TRUE : PETSC_FALSE;
+  element->property = iga->property ? PETSC_TRUE : PETSC_FALSE;
   { /* */
-    element->rational = iga->rational ? PETSC_TRUE : PETSC_FALSE;
-    element->geometry = iga->geometry ? PETSC_TRUE : PETSC_FALSE;
-    element->property = iga->property ? PETSC_TRUE : PETSC_FALSE;
+    PetscInt nen = element->nen;
+    PetscInt nsd = iga->geometry ? iga->geometry : iga->dim;
+    PetscInt npd = iga->property;
+    if (element->nsd != nsd) {
+      element->nsd = nsd;
+      ierr = PetscFree(element->geometryX);CHKERRQ(ierr);
+      ierr = PetscMalloc1(nen*nsd,&element->geometryX);CHKERRQ(ierr);
+    }
+    if (element->npd != npd) {
+      element->npd = npd;
+      ierr = PetscFree(element->propertyA);CHKERRQ(ierr);
+      ierr = PetscMalloc1(nen*npd,&element->propertyA);CHKERRQ(ierr);
+    }
+    ierr = PetscMemzero(element->rationalW,sizeof(PetscReal)*nen);CHKERRQ(ierr);
+    ierr = PetscMemzero(element->geometryX,sizeof(PetscReal)*nen*nsd);CHKERRQ(ierr);
+    ierr = PetscMemzero(element->propertyA,sizeof(PetscReal)*nen*npd);CHKERRQ(ierr);
   }
   { /* */
     PetscInt q,i;
     point->dim = element->dim;
     point->nsd = element->nsd;
     point->npd = element->npd;
+    point->rational = element->rational ? element->rationalW : NULL;
+    point->geometry = element->geometry ? element->geometryX : NULL;
+    point->property = element->property ? element->propertyA : NULL;
   }
   PetscFunctionReturn(0);
 }
 PetscBool IGAElementNextPoint(IGAElement element,IGAPoint point)
 {
   PetscInt nen = point->nen;
-  PetscInt nsd = point->nsd;
   PetscInt dim = point->dim;
   PetscInt dim2 = dim*dim;
   PetscInt dim3 = dim*dim2;
 
  start:
 
-  point->rational = element->rationalW;
-  point->geometry = element->geometryX;
-  point->property = element->propertyA;
-  if (!element->rational) point->rational = NULL;
-  if (!element->geometry) point->geometry = NULL;
-  if (!element->property) point->property = NULL;
-
   point->point    = element->point;
   point->weight   = element->weight;
   point->detJac   = element->detJac;
   point->detS     = element->detS;
   point->normal   = element->normal;
 
-  if (element->geometry && dim == nsd) { /* XXX */
+  if (element->geometry && dim == point->nsd) { /* XXX */
     point->shape[0] = element->shape[0];
     point->shape[1] = element->shape[1];
     point->shape[2] = element->shape[2];
 @*/
 PetscErrorCode IGASetGeometryDim(IGA iga,PetscInt dim)
 {
+  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidLogicalCollectiveInt(iga,dim,2);
     SETERRQ2(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
              "Cannot change IGA geometry dim to %D after it was set to %D",dim,iga->geometry);
   if (iga->geometry == dim) PetscFunctionReturn(0);
+  ierr = PetscFree(iga->geometryX);CHKERRQ(ierr);
+  ierr = PetscFree(iga->rationalW);CHKERRQ(ierr);
+  iga->rational = PETSC_FALSE;
   iga->geometry = dim;
   iga->setup = PETSC_FALSE;
   PetscFunctionReturn(0);
 @*/
 PetscErrorCode IGASetPropertyDim(IGA iga,PetscInt dim)
 {
+  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidLogicalCollectiveInt(iga,dim,2);
     SETERRQ2(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
              "Cannot change IGA property dim to %D after it was set to %D",dim,iga->property);
   if (iga->property == dim) PetscFunctionReturn(0);
+  ierr = PetscFree(iga->propertyA);CHKERRQ(ierr);
   iga->property = dim;
   iga->setup = PETSC_FALSE;
   PetscFunctionReturn(0);

src/petigapoint.c

   PetscValidPointer(point,1);
   point->count =  0;
   point->index = -1;
+  point->rational = NULL;
+  point->geometry = NULL;
+  point->property = NULL;
   ierr = IGAPointFreeWork(point);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
   point->dim = element->dim;
   point->nsd = element->nsd;
   point->npd = element->npd;
+  point->rational = element->rational ? element->rationalW : NULL;
+  point->geometry = element->geometry ? element->geometryX : NULL;
+  point->property = element->property ? element->propertyA : NULL;
   { /* */
     size_t MAX_WORK_VEC = sizeof(point->wvec)/sizeof(PetscScalar*);
     size_t MAX_WORK_MAT = sizeof(point->wmat)/sizeof(PetscScalar*);
-    size_t i, n = point->nen * point->dof;
+    size_t i, nv = element->nen * element->dof, nm = nv*nv;
     for (i=0; i<MAX_WORK_VEC; i++)
-      {ierr = PetscMalloc1(n,&point->wvec[i]);CHKERRQ(ierr);}
+      {ierr = PetscMalloc1(nv,&point->wvec[i]);CHKERRQ(ierr);}
     for (i=0; i<MAX_WORK_MAT; i++)
-      {ierr = PetscMalloc1(n*n,&point->wmat[i]);CHKERRQ(ierr);}
+      {ierr = PetscMalloc1(nm,&point->wmat[i]);CHKERRQ(ierr);}
   }
   PetscFunctionReturn(0);
 }

test/GeometryMap.c

   ierr = IGASetOrder(iga,3);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
+  ierr = IGASetGeometryDim(iga,dim);CHKERRQ(ierr);
   iga->rational = PETSC_TRUE;
-  ierr = IGASetGeometryDim(iga,dim);CHKERRQ(ierr);
   ierr = PetscMalloc(3*3*(dim-1)*sizeof(PetscReal),&iga->rationalW);CHKERRQ(ierr);
   ierr = PetscMalloc(3*3*(dim-1)*dim*sizeof(PetscReal),&iga->geometryX);CHKERRQ(ierr);
   {