1. Lisandro Dalcin
  2. PetIGA

Commits

Lisandro Dalcin  committed 82de52a

Enhancements to boundary integrals implementation

  • Participants
  • Parent commits a686184
  • Branches default

Comments (0)

Files changed (3)

File include/petiga.h

View file
  • Ignore whitespace
   PetscReal *weight;  /* [nqp]                */
   PetscReal *point;   /* [nel][nqp]           */
   PetscReal *value;   /* [nel][nqp][nen][d+1] */
+
+  PetscReal  bnd_detJ[2];
+  PetscReal  bnd_weight[2];
+  PetscReal  bnd_point[2];
+  PetscReal *bnd_value[2]; /* [nen][d+1] */
 };
 PETSC_EXTERN PetscErrorCode IGABasisCreate(IGABasis *basis);
 PETSC_EXTERN PetscErrorCode IGABasisDestroy(IGABasis *basis);

File src/petigabasis.c

View file
  • Ignore whitespace
   ierr = PetscFree(basis->weight);CHKERRQ(ierr);
   ierr = PetscFree(basis->point);CHKERRQ(ierr);
   ierr = PetscFree(basis->value);CHKERRQ(ierr);
+  ierr = PetscFree(basis->bnd_value[0]);CHKERRQ(ierr);
+  ierr = PetscFree(basis->bnd_value[1]);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   basis->point  = point;
   basis->value  = value;
 
+  {
+    PetscInt  k0 = span[0], k1 = span[nel-1];
+    PetscReal u0 = U[k0],   u1 = U[k1+1];
+    ierr = PetscMalloc1(nen*ndr,PetscReal,&basis->bnd_value[0]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen*ndr,PetscReal,&basis->bnd_value[1]);CHKERRQ(ierr);
+    basis->bnd_detJ  [0] = 1.0; basis->bnd_detJ  [1] = 1.0;
+    basis->bnd_weight[0] = 1.0; basis->bnd_weight[1] = 1.0;
+    basis->bnd_point [0] =  u0; basis->bnd_point [1] = u1;
+    IGA_DersBasisFuns(k0,u0,p,d,U,basis->bnd_value[0]);
+    IGA_DersBasisFuns(k1,u1,p,d,U,basis->bnd_value[1]);
+  }
+
   PetscFunctionReturn(0);
 }
 

File src/petigaelem.c

View file
  • Ignore whitespace
   PetscFunctionReturn(0);
 }
 
+#define IGA_Quadrature_BNDS(ID,BD,i,s) \
+  1,&BD[i]->bnd_point[s],&BD[i]->bnd_weight[s],&BD[i]->bnd_detJ[s]
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementBuildBoundaryQuadrature"
 PetscErrorCode IGAElementBuildBoundaryQuadrature(IGAElement element,PetscInt dir,PetscInt side)
   {
     IGABasis *BD = element->BD;
     PetscInt *ID = element->ID;
-
-    PetscReal U[2],one=1;
-    IGAAxisGetLimits(element->parent->axis[dir],&U[0],&U[1]);
-
     switch (element->dim) {
     case 3:
       switch (dir) {
-      case 0: IGA_Quadrature_3D(1,&U[side],&one,&one,
+      case 0: IGA_Quadrature_3D(IGA_Quadrature_BNDS(ID,BD,0,side),
                                 IGA_Quadrature_ARGS(ID,BD,1),
                                 IGA_Quadrature_ARGS(ID,BD,2),
                                 element->weight,element->detJac,
                                 element->point,element->scale); break;
       case 1: IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
-                                1,&U[side],&one,&one,
+                                IGA_Quadrature_BNDS(ID,BD,1,side),
                                 IGA_Quadrature_ARGS(ID,BD,2),
                                 element->weight,element->detJac,
                                 element->point,element->scale); break;
       case 2: IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
                                 IGA_Quadrature_ARGS(ID,BD,1),
-                                1,&U[side],&one,&one,
+                                IGA_Quadrature_BNDS(ID,BD,2,side),
                                 element->weight,element->detJac,
                                 element->point,element->scale); break;
       }
       break;
     case 2:
       switch (dir) {
-      case 0: IGA_Quadrature_2D(1,&U[side],&one,&one,
+      case 0: IGA_Quadrature_2D(IGA_Quadrature_BNDS(ID,BD,0,side),
                                 IGA_Quadrature_ARGS(ID,BD,1),
                                 element->weight,element->detJac,
                                 element->point,element->scale); break;
       case 1: IGA_Quadrature_2D(IGA_Quadrature_ARGS(ID,BD,0),
-                                1,&U[side],&one,&one,
+                                IGA_Quadrature_BNDS(ID,BD,1,side),
                                 element->weight,element->detJac,
                                 element->point,element->scale); break;
       }
       break;
     case 1:
       switch (dir) {
-      case 0: IGA_Quadrature_1D(1,&U[side],&one,&one,
+      case 0: IGA_Quadrature_1D(IGA_Quadrature_BNDS(ID,BD,0,side),
                                 element->weight,element->detJac,
                                 element->point,element->scale); break;
       }
   PetscFunctionReturn(0);
 }
 
-EXTERN_C_BEGIN
-extern void IGA_DersBasisFuns(PetscInt i,PetscReal u,PetscInt p,PetscInt d,const PetscReal U[],PetscReal N[]);
-EXTERN_C_END
+#define IGA_BasisFuns_BNDS(ID,BD,i,s) \
+  1,BD[i]->nen,BD[i]->d,BD[i]->bnd_value[s]
 
 EXTERN_C_BEGIN
 extern void IGA_GetNormal(PetscInt dim,PetscInt dir,PetscInt side,const PetscReal F[],PetscReal *dS,PetscReal N[]);
     IGABasis *BD  = element->BD;
     PetscInt *ID  = element->ID;
     PetscReal **N = element->basis;
-
-    PetscInt k = element->parent->axis[dir]->span[element->ID[dir]];
-    PetscInt p = element->parent->axis[dir]->p;
-    PetscInt d = element->parent->order;
-    PetscReal *U = element->parent->axis[dir]->U;
-    PetscInt nen = BD[dir]->nen;
-    PetscReal value[nen*(d+1)];
-
-    PetscReal u[2];
-    IGAAxisGetLimits(element->parent->axis[dir],&u[0],&u[1]);
-    IGA_DersBasisFuns(k,u[side],p,d,U,value);
-
     switch (element->dim) {
     case 3:
       switch (dir) {
       case 0: IGA_BasisFuns_3D(order,
                                element->rational, element->rationalW,
-                               1,BD[0]->nen,BD[0]->d,value,
+                               IGA_BasisFuns_BNDS(ID,BD,0,side),
                                IGA_BasisFuns_ARGS(ID,BD,1),
                                IGA_BasisFuns_ARGS(ID,BD,2),
                                N[0],N[1],N[2],N[3]); break;
       case 1: IGA_BasisFuns_3D(order,
                                element->rational,element->rationalW,
                                IGA_BasisFuns_ARGS(ID,BD,0),
-                               1,BD[1]->nen,BD[1]->d,value,
+                               IGA_BasisFuns_BNDS(ID,BD,1,side),
                                IGA_BasisFuns_ARGS(ID,BD,2),
                                N[0],N[1],N[2],N[3]); break;
       case 2: IGA_BasisFuns_3D(order,
                                element->rational,element->rationalW,
                                IGA_BasisFuns_ARGS(ID,BD,0),
                                IGA_BasisFuns_ARGS(ID,BD,1),
-                               1,BD[2]->nen,BD[2]->d,value,
+                               IGA_BasisFuns_BNDS(ID,BD,2,side),
                                N[0],N[1],N[2],N[3]); break;
       }
       break;
       switch (dir) {
       case 0: IGA_BasisFuns_2D(order,
                                element->rational,element->rationalW,
-                               1,BD[0]->nen,BD[0]->d,value,
+                               IGA_BasisFuns_BNDS(ID,BD,0,side),
                                IGA_BasisFuns_ARGS(ID,BD,1),
                                N[0],N[1],N[2],N[3]); break;
       case 1: IGA_BasisFuns_2D(order,
                                element->rational,element->rationalW,
                                IGA_BasisFuns_ARGS(ID,BD,0),
-                               1,BD[1]->nen,BD[1]->d,value,
+                               IGA_BasisFuns_BNDS(ID,BD,1,side),
                                N[0],N[1],N[2],N[3]); break;
       }
       break;
       switch (dir) {
       case 0: IGA_BasisFuns_1D(order,
                                element->rational,element->rationalW,
-                               1,BD[0]->nen,BD[0]->d,value,
+                               IGA_BasisFuns_BNDS(ID,BD,0,side),
                                N[0],N[1],N[2],N[3]); break;
       }
       break;
   if (!element->geometry) {
     PetscInt q, nqp = element->nqp;
     PetscInt i, dim = element->dim;
-    PetscReal *n  = element->normal;
+    PetscReal *n = element->normal;
     for (q=0; q<nqp; q++)
       for (i=0; i<dim; i++)
-        n[q*dim+i] = (dir==i) ? ((side==0)?-1.0:1.0) : 0.0;
+        n[q*dim+i] = (dir==i) ? ((side==0) ? -1.0 : 1.0) : 0.0;
   }
   PetscFunctionReturn(0);
 }