Commits

Lisandro Dalcin  committed 1b5cd15

Changes to the handling of collocation

  • Participants
  • Parent commits e1c84e0

Comments (0)

Files changed (4)

File include/petiga.h

   IGARule     rule[3];
   IGABoundary boundary[3][2];
 
-  IGABasis    elem_basis[3];
-  IGAElement  elem_iterator;
-
-  /* stuff added for collocation */
-  PetscBool   collocation;
-  IGABasis    node_basis[3];
-  IGAElement  node_iterator;
-
+  IGABasis    basis[3];
+  IGAElement  iterator;
 
   PetscInt    geometry;
   PetscBool   rational;
   Vec        natural;
   VecScatter n2g,g2n;
 
+  /* stuff added for collocation */
+  PetscBool   collocation;
+  IGABasis    node_basis[3];
+  IGAElement  node_iterator;
 };
 
 PETSC_EXTERN PetscClassId IGA_CLASSID;

File src/petiga.c

   for (i=0; i<3; i++) {
     ierr = IGAAxisCreate(&iga->axis[i]);CHKERRQ(ierr);
     ierr = IGARuleCreate(&iga->rule[i]);CHKERRQ(ierr);
-    ierr = IGABasisCreate(&iga->elem_basis[i]);CHKERRQ(ierr);
-    ierr = IGABasisCreate(&iga->node_basis[i]);CHKERRQ(ierr);
+    ierr = IGABasisCreate(&iga->basis[i]);CHKERRQ(ierr);
     ierr = IGABoundaryCreate(&iga->boundary[i][0]);CHKERRQ(ierr);
     ierr = IGABoundaryCreate(&iga->boundary[i][1]);CHKERRQ(ierr);
     iga->proc_sizes[i] = -1;
   }
-  ierr = IGAElementCreate(&iga->elem_iterator);CHKERRQ(ierr);
-  ierr = IGAElementCreate(&iga->node_iterator);CHKERRQ(ierr);
+  ierr = IGAElementCreate(&iga->iterator);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
   for (i=0; i<3; i++) {
     ierr = IGAAxisDestroy(&iga->axis[i]);CHKERRQ(ierr);
     ierr = IGARuleDestroy(&iga->rule[i]);CHKERRQ(ierr);
-    ierr = IGABasisDestroy(&iga->elem_basis[i]);CHKERRQ(ierr);
-    ierr = IGABasisDestroy(&iga->node_basis[i]);CHKERRQ(ierr);
+    ierr = IGABasisDestroy(&iga->basis[i]);CHKERRQ(ierr);
     ierr = IGABoundaryDestroy(&iga->boundary[i][0]);CHKERRQ(ierr);
     ierr = IGABoundaryDestroy(&iga->boundary[i][1]);CHKERRQ(ierr);
   }
-  ierr = IGAElementDestroy(&iga->elem_iterator);CHKERRQ(ierr);
+  ierr = IGAElementDestroy(&iga->iterator);CHKERRQ(ierr);
+
+  for (i=0; i<3; i++) /* collocation */
+    {ierr = IGABasisDestroy(&iga->node_basis[i]);CHKERRQ(ierr);}
   ierr = IGAElementDestroy(&iga->node_iterator);CHKERRQ(ierr);
 
   ierr = IGAReset(iga);CHKERRQ(ierr);
     {ierr = VecDestroy(&iga->vwork[--iga->nwork]);CHKERRQ(ierr);}
   ierr = DMDestroy(&iga->node_dm);CHKERRQ(ierr);
 
-  ierr = IGAElementReset(iga->elem_iterator);CHKERRQ(ierr);
+  ierr = IGAElementReset(iga->iterator);CHKERRQ(ierr);
   ierr = IGAElementReset(iga->node_iterator);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidLogicalCollectiveBool(iga,collocation,2);
-  if (!iga->collocation && collocation) {
+  if (collocation && !iga->collocation) {
     PetscMPIInt size = 1;
     PetscInt i, dim = (iga->dim > 0) ? iga->dim : 3;
     PetscBool periodic = PETSC_FALSE;
     if (periodic) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_SUP,
 			  "Collocation not supported with periodicity");
   }
+  if (collocation && iga->setup) { /* collocation */
+    PetscInt i;
+    for (i=0; i<3; i++) {
+      ierr = IGABasisDestroy(&iga->node_basis[i]);CHKERRQ(ierr);
+      ierr = IGABasisCreate(&iga->node_basis[i]);CHKERRQ(ierr);
+      ierr = IGABasisInitCollocation(iga->node_basis[i],iga->axis[i],iga->order);CHKERRQ(ierr);
+    }
+    ierr = IGAElementDestroy(&iga->node_iterator);CHKERRQ(ierr);
+    ierr = IGAElementCreate(&iga->node_iterator);CHKERRQ(ierr);
+    iga->node_iterator->collocation = PETSC_TRUE;
+    ierr = IGAElementInit(iga->node_iterator,iga);CHKERRQ(ierr);
+  }
   iga->collocation = collocation;
   PetscFunctionReturn(0);
 }
   IGACheckSetUp(iga,1);
   if (i < 0)         SETERRQ1(((PetscObject)iga)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Index %D must be nonnegative",i);
   if (i >= iga->dim) SETERRQ2(((PetscObject)iga)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Index %D, but dimension %D",i,iga->dim);
-  *basis = iga->elem_basis[i];
+  *basis = iga->basis[i];
   PetscFunctionReturn(0);
 }
 
   for (i=0; i<3; i++) {
     if (i >= iga->dim) {ierr = IGARuleReset(iga->rule[i]);CHKERRQ(ierr);}
     if (iga->rule[i]->nqp < 1) {ierr = IGARuleInit(iga->rule[i],iga->axis[i]->p + 1);CHKERRQ(ierr);}
-    ierr = IGABasisInitQuadrature(iga->elem_basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);
+    ierr = IGABasisInitQuadrature(iga->basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);
   }
-  iga->elem_iterator->collocation = PETSC_FALSE;
-  ierr = IGAElementInit(iga->elem_iterator,iga);CHKERRQ(ierr);
-
-  for (i=0; i<3; i++) {
-    ierr = IGABasisInitCollocation(iga->node_basis[i],iga->axis[i],iga->order);CHKERRQ(ierr);
-  }
-  iga->node_iterator->collocation = PETSC_TRUE;
-  ierr = IGAElementInit(iga->node_iterator,iga);CHKERRQ(ierr);
+  ierr = IGAElementInit(iga->iterator,iga);CHKERRQ(ierr);
 
   ierr = IGASetUp_View(iga);CHKERRQ(ierr);
+
+  if (iga->collocation) {ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr);}  /* collocation */
+
   PetscFunctionReturn(0);
 }
 

File src/petigaelem.c

   element->dim = iga->dim;
   element->nsd = iga->geometry ? iga->geometry : iga->dim;
   element->npd = iga->property ? iga->property : 0;
-  if (!element->collocation) {
-    start = iga->elem_start;
-    width = iga->elem_width;
-    sizes = iga->elem_sizes;
-    element->BD = iga->elem_basis;
-  } else {
+
+  start = iga->elem_start;
+  width = iga->elem_width;
+  sizes = iga->elem_sizes;
+  element->BD = iga->basis;
+  if (PetscUnlikely(element->collocation)) { /* collocation */
     start = iga->node_lstart;
     width = iga->node_lwidth;
     sizes = iga->node_sizes;
 
     ierr = PetscMemzero(element->normal,sizeof(PetscReal)*nqp*dim);CHKERRQ(ierr);
   }
-  if (!element->collocation) {
-    element->neq    = element->nen;
-    element->rowmap = element->mapping;
-    element->colmap = element->mapping;
-  } else {
+  element->neq    = element->nen;
+  element->rowmap = element->mapping;
+  element->colmap = element->mapping;
+  if (PetscUnlikely(element->collocation)) { /* collocation */
     element->neq    = 1;
     element->rowmap = &element->index;
     element->colmap = element->mapping;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidPointer(element,2);
-  if (!iga->collocation)
-    *element = iga->elem_iterator;
-  else
+  *element = iga->iterator;
+  if (PetscUnlikely(iga->collocation)) /* collocation */
     *element = iga->node_iterator;
   PetscFunctionReturn(0);
 }
   PetscValidPointer(element,1);
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
+  if (PetscUnlikely(element->collocation)) goto collocation;
   element->nfix  = 0;
   element->nflux = 0;
-  if (!element->collocation) {
+  {
     IGAAxis  *AX = element->parent->axis;
     PetscInt *ID = element->ID;
     PetscInt i,dim = element->dim;
       if (ID[i] == 0 && !w) BuildFix(element,i,0);
       if (ID[i] == e && !w) BuildFix(element,i,1);
     }
-  } else {
+  } 
+  PetscFunctionReturn(0);
+ collocation: 
+  element->nfix  = 0;
+  element->nflux = 0;
+  {
     PetscInt A0[3] = {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT};
     PetscInt A1[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
     {
   PetscValidPointer(element,1);
   PetscValidScalarPointer(K,2);
   PetscValidScalarPointer(F,3);
-  if (!element->collocation) {
+  if (PetscUnlikely(element->collocation)) goto collocation;
+  {
     PetscInt M = element->neq * element->dof;
     PetscInt N = element->nen * element->dof;
     PetscInt f,n;
       K[k*N+k] = 1.0;
       F[k]     = v;
     }
-  } else {
+  }
+  PetscFunctionReturn(0);
+ collocation:
+  {
     PetscInt M = element->neq * element->dof;
     PetscInt N = element->nen * element->dof;
     PetscInt f,n;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   PetscValidScalarPointer(F,2);
-  if (!element->collocation) {
+  if (PetscUnlikely(element->collocation)) goto collocation;
+  {
     PetscInt f,n;
     n = element->nflux;
     for (f=0; f<n; f++) {
       PetscScalar u = element->ufix[f];
       F[k] = u - v;
     }
-  } else {
+  }
+  PetscFunctionReturn(0);
+ collocation: 
+  {
     PetscInt f,n;
     n = element->nfix;
     for (f=0; f<n; f++) {
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   PetscValidScalarPointer(J,2);
-  if (!element->collocation) {
+  if (PetscUnlikely(element->collocation)) goto collocation;
+  {
     PetscInt M = element->neq * element->dof;
     PetscInt N = element->nen * element->dof;
     PetscInt f,n;
       for (j=0; j<N; j++) J[k*N+j] = 0.0;
       J[k*N+k] = 1.0;
     }
-  } else {
+  }
+  PetscFunctionReturn(0);
+ collocation:
+  {
     PetscInt M = element->neq * element->dof;
     PetscInt N = element->nen * element->dof;
     PetscInt f,n;

File src/petigamat.c

     *first = k - s - n;
   }
 
-  if (iga->collocation) {
+  if (PetscUnlikely(iga->collocation)) { /* collocation */
     PetscInt offset = iga->node_basis[dir]->offset[i];
     *first = offset;
     *last  = offset + p;