Commits

Lisandro Dalcin committed e2df752

Update demos to use typeof for VLA casts

  • Participants
  • Parent commits 38c7bce

Comments (0)

Files changed (11)

File demo/AdvectionDiffusion.c

   PetscReal adv[3];
 } AppCtx;
 
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   AppCtx *user = (AppCtx *)ctx;
 
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
-
-  const PetscReal *N0,(*N1)[dim];
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
-  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+  const PetscReal *N0        = (typeof(N0)) p->shape[0];
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
 
   PetscInt a,b,i;
   for (a=0; a<nen; a++) {
     for (b=0; b<nen; b++) {
-      PetscScalar diffusion = 0.0;
+      PetscScalar diffusion = DOT(dim,N1[a],N1[b]);
       PetscScalar advection = 0.0;
-      for (i=0; i<dim; i++){
-        diffusion += N1[a][i]*N1[b][i];
-	advection += user->adv[i]*N1[b][i];
-      }
+      for (i=0; i<dim; i++)
+        advection += user->adv[i]*N1[b][i];
       advection *= N0[a];
-      K[a*nen+b] = diffusion + advection; 
+      K[a*nen+b] = diffusion + advection;
     }
     F[a] = 0.0;
   }
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
-  
+
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   IGABoundary bnd;
   for (i=0; i<dim; i++) {

File demo/BoundaryIntegral.c

 {
   PetscInt dim = p->dim;
   PetscInt nen = p->nen;
-  const PetscReal (*N1)[dim]; //*((PetscReal**)&N1) = p->shape[1];
-  IGAPointGetShapeFuns(p,1,(const PetscReal **)&N1);
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
   PetscInt a,b;
   for (a=0; a<nen; a++) {
     for (b=0; b<nen; b++)
 PetscErrorCode Neumann(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt nen = p->nen;
-  const PetscReal *N0 = p->shape[0];
+  const PetscReal *N0 = (typeof(N0)) p->shape[0];
   PetscInt a;
   for (a=0; a<nen; a++)
     F[a] = N0[a] * 1.0;
 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
-
-  const PetscReal (*N2)[dim][dim]; //*((PetscReal**)&N2) = p->shape[2];
-  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+  const PetscReal (*N2)[dim][dim] = (typeof(N2)) p->shape[2];
 
   PetscInt a;
   for (a=0; a<nen; a++)
 PetscErrorCode DirichletCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt nen = p->nen;
-  const PetscReal *N0 = p->shape[0];
+  const PetscReal *N0 = (typeof(N0)) p->shape[0];
+
   PetscInt a;
   for (a=0; a<nen; a++)
     K[a] = N0[a];
 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
-  const PetscReal (*N1)[dim];
-  const PetscReal *normal = p->normal;
-  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
+  const PetscReal *normal    = p->normal;
   PetscInt a;
   for (a=0; a<nen; a++)
     K[a] = DOT(dim,N1[a],normal);

File demo/Bratu.c

-/* 
+/*
    This code solves the steady and unsteady Bratu equation. It also
    demonstrates how the user-specified routines, here the Function and
    Jacobian routines, can be implemented in Fortran (see BratuFJ.F90)
 */
 #include "petiga.h"
 
-typedef struct { 
+typedef struct {
   PetscReal lambda;
 } AppCtx;
 
 
   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
-  
-  
+
+
   PetscBool steady = PETSC_TRUE;
   PetscReal lambda = 6.80;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Bratu Options","IGA");CHKERRQ(ierr);
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
-
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-  IGABoundary bnd;
   PetscInt dir,side;
   for (dir=0; dir<3; dir++) {
     for (side=0; side<2; side++) {
+      IGABoundary bnd;
+      PetscInt    field = 0;
+      PetscScalar value = 0.0;
       ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,field,value);CHKERRQ(ierr);
     }
   }
-
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-
   PetscInt dim;
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   if (dim < 1) {ierr = IGASetDim(iga,dim=2);CHKERRQ(ierr);}

File demo/L2Projection.c

 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen;
-  IGAPointGetSizes(p,0,&nen,0);
+  PetscInt nen = p->nen;
 
   PetscReal x[3] = {0,0,0};
   IGAPointFormPoint(p,x);
   PetscScalar f = Function(x[0],x[1],x[2]);
 
-  const PetscReal *N;
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N);
+  const PetscReal *N = (typeof(N)) p->shape[0];
 
   PetscInt a,b;
   for (a=0; a<nen; a++) {

File demo/Laplace.c

 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
-
-  const PetscReal (*N1)[dim];
-  IGAPointGetShapeFuns(p,1,(const PetscReal **)&N1);
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
 
   PetscInt a,b;
   for (a=0; a<nen; a++) {
 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
-
-  const PetscReal (*N2)[dim][dim];
-  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+  const PetscReal (*N2)[dim][dim] = (typeof(N2)) p->shape[2];
 
   PetscInt a;
   for (a=0; a<nen; a++)

File demo/LoggChallenge.c

 {
   PetscReal x = p->point[0];
   PetscReal y = p->point[1];
+  PetscReal f = Forcing(x,y);
   PetscReal *N0 = p->shape[0];
-  PetscReal (*N1)[2] = (PetscReal (*)[2]) p->shape[1];
+  PetscReal (*N1)[2] = (PetscReal(*)[2]) p->shape[1];
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na   = N0[a];
       PetscReal Nb_y = N1[b][1];
       K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y;
     }
-    F[a] = Na * Forcing(x,y);
+    F[a] = Na * f;
   }
   return 0;
 }

File demo/Neumann.c

 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
+  const PetscReal *N0        = (typeof(N0)) p->shape[0];
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
 
   PetscReal x[3] = {0,0,0};
   IGAPointFormPoint(p,x);
   PetscReal f = Forcing(x);
 
-  const PetscReal *N0,(*N1)[dim];
-  IGAPointGetShapeFuns(p,0,(const PetscReal **)&N0);
-  IGAPointGetShapeFuns(p,1,(const PetscReal **)&N1);
-
   PetscInt a,b;
   for (a=0; a<nen; a++) {
     for (b=0; b<nen; b++)
 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
+  const PetscReal (*N2)[dim][dim] = (typeof(N2)) p->shape[2];
 
   PetscReal x[3] = {0,0,0};
   IGAPointFormPoint(p,x);
   PetscReal f = Forcing(x);
 
-  const PetscReal (*N2)[dim][dim];
-  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
-
   PetscInt a;
   for (a=0; a<nen; a++)
     K[a] += -DEL2(dim,N2[a]);

File demo/PatternFormation.c

 /*
   This code solves the following system of PDEs
-  
+
      u_t = D_1 \nabla^2 u + f(u,v)
      v_t = D_2 \nabla^2 v + g(u,v)
 
   PetscReal tau1  = user->tau1;
   PetscReal tau2  = user->tau2;
 
-  PetscInt nen;
-  IGAPointGetSizes(p,0,&nen,0);
+  PetscInt nen = p->nen;
   PetscScalar (*R)[2] = (PetscScalar (*)[2])Re;
 
   PetscScalar uv_t[2],uv_0[2],uv_1[2][2];
   PetscReal f = alpha*u*(1-tau1*v*v) + v*(1-tau2*u);
   PetscReal g = beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);
 
-  const PetscReal *N0,(*N1)[2];
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
-  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+  const PetscReal *N0      = p->shape[0];
+  const PetscReal (*N1)[2] = (const PetscReal(*)[2]) p->shape[1];
 
-  PetscInt  a;
+  PetscInt a;
   for (a=0; a<nen; a++) {
     PetscReal Na   = N0[a];
     PetscReal Na_x = N1[a][0];
   PetscReal tau1  = user->tau1;
   PetscReal tau2  = user->tau2;
 
-  PetscInt nen;
-  IGAPointGetSizes(p,0,&nen,0);
-  PetscScalar (*K)[2][nen][2] = (PetscScalar (*)[2][nen][2])Ke;
+  PetscInt nen = p->nen;
+  PetscScalar (*K)[2][nen][2] = (typeof(K)) Ke;
 
   PetscReal f_u=0,f_v=0;
   PetscReal g_u=0,g_v=0;

File demo/Poisson.c

 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
-
-  const PetscReal *N0,(*N1)[dim];
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
-  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+  const PetscReal *N0        = (typeof(N0)) p->shape[0];
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
 
   PetscInt a,b;
   for (a=0; a<nen; a++) {
 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
-
-  const PetscReal (*N2)[dim][dim];
-  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+  const PetscReal (*N2)[dim][dim] = (typeof(N2)) p->shape[2];
 
   PetscInt a;
   for (a=0; a<nen; a++)

File demo/Richards.c

 #define SQ(A) ((A)*(A))
 
 typedef struct {
-  PetscScalar m,lambda,kappa,NGr; 
+  PetscScalar m,lambda,kappa,NGr;
   PetscScalar S0,Sin,d0,penalty;
   PetscBool   phase_field;
 } AppCtx;
 PetscErrorCode L2Projection(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   AppCtx *user = (AppCtx *)ctx;
-  PetscInt a,b,dim,nen;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,NULL,NULL);
+
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
 
   PetscReal x[dim];
   IGAPointFormPoint(p,x);
-  
+
   PetscScalar f = user->S0;
   PetscScalar d = x[dim-1];
   if(d < user->d0) f += SQ(1.-SQ(d/user->d0))*(user->Sin-user->S0);
 
-  const PetscReal *N;
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N);
+  const PetscReal *N = (typeof(N)) p->shape[0];
+
+  PetscInt a,b;
   for(a=0; a<nen; a++){
-    for(b=0; b<nen; b++)  K[b*nen+a] = N[a]*N[b];
+    for(b=0; b<nen; b++) K[b*nen+a] = N[a]*N[b];
     F[a] = N[a]*f;
   }
   return 0;
 #undef  __FUNCT__
 #define __FUNCT__ "Residual"
 PetscErrorCode Residual(IGAPoint p,PetscReal dt,
-			   PetscReal shift,const PetscScalar *V,
-			   PetscReal t,const PetscScalar *U,
-			   PetscScalar *R,void *ctx)
+                           PetscReal shift,const PetscScalar *V,
+                           PetscReal t,const PetscScalar *U,
+                           PetscScalar *R,void *ctx)
 {
   AppCtx *user = (AppCtx *)ctx;
 
-  PetscInt a,i,nen,dim;
-  IGAPointGetSizes(p,&nen,0,0);
-  IGAPointGetDims(p,&dim,NULL,NULL);
+  PetscInt a,i;
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
 
   PetscScalar S_t,S,delS=0;
   PetscScalar S1[dim],S2[dim][dim];
   for(i=0;i<dim;i++) delS += S2[i][i];
 
   S = PetscMax(1000*PETSC_MACHINE_EPSILON,S);
-  
-  const PetscReal *N0,(*N1)[dim],(*N2)[dim][dim];
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
-  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
-  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+
+  const PetscReal  *N0            = (typeof(N0)) p->shape[0];
+  const PetscReal (*N1)[dim]      = (typeof(N1)) p->shape[1];
+  const PetscReal (*N2)[dim][dim] = (typeof(N2)) p->shape[2];
 
   PetscScalar kappa     = user->kappa;
   PetscScalar lambda    = user->lambda;
   CapillaryPressure(S,lambda,kappa,NULL,&dJdS);
 
   delS *= rNGr*rNGr*rNGr; // include the NGr^-3 into delS
+
   for (a=0; a<nen; a++) {
     PetscScalar gN_dot_gS=0,delN=0,gkD_dot_gN=0;
     for(i=0;i<dim;i++) {
 #undef  __FUNCT__
 #define __FUNCT__ "Jacobian"
 PetscErrorCode Jacobian(IGAPoint p,PetscReal dt,
-			PetscReal shift,const PetscScalar *V,
-			PetscReal t,const PetscScalar *U,
-			PetscScalar *J,void *ctx)
+                        PetscReal shift,const PetscScalar *V,
+                        PetscReal t,const PetscScalar *U,
+                        PetscScalar *J,void *ctx)
 {
-  
+
   return 0;
 }
 
 #undef  __FUNCT__
 #define __FUNCT__ "ZeroFluxResidual"
 PetscErrorCode ZeroFluxResidual(IGAPoint p,PetscReal dt,
-				PetscReal shift,const PetscScalar *V,
-				PetscReal t,const PetscScalar *U,
-				PetscScalar *R,void *ctx)
+                                PetscReal shift,const PetscScalar *V,
+                                PetscReal t,const PetscScalar *U,
+                                PetscScalar *R,void *ctx)
 {
   AppCtx *user = (AppCtx *)ctx;
 
-  PetscInt a,i,nen,dim;
-  IGAPointGetSizes(p,&nen,0,0);
-  IGAPointGetDims(p,&dim,NULL,NULL);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
 
-  PetscScalar *n = p->normal;
-  
   PetscScalar S1[dim];
   IGAPointFormGrad(p,U,&S1[0]);
 
-  const PetscReal *N0;
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+  const PetscReal *N0 = (typeof(N0)) p->shape[0];
+  const PetscReal *n = p->normal;
 
+  PetscInt a,i;
   for (a=0; a<nen; a++) {
-    PetscScalar n_gS = 0.;
+    PetscScalar n_gS = 0.0;
     for(i=0;i<dim;i++) n_gS += n[i]*S1[i];
     R[a] = user->penalty*N0[a]*n_gS;
   }
 #undef  __FUNCT__
 #define __FUNCT__ "ZeroFluxJacobian"
 PetscErrorCode ZeroFluxJacobian(IGAPoint p,PetscReal dt,
-				PetscReal shift,const PetscScalar *V,
-				PetscReal t,const PetscScalar *U,
-				PetscScalar *J,void *ctx)
+                                PetscReal shift,const PetscScalar *V,
+                                PetscReal t,const PetscScalar *U,
+                                PetscScalar *J,void *ctx)
 {
 
   return 0;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
   AppCtx user;
-  user.m       = 7; 
+  user.m       = 7;
   user.kappa   = 50;
   user.lambda  = 4;
   user.NGr     = 50;
   }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
-  
+
   // Boundary conditions
   for(dir=0;dir<dim;dir++)
     for(side=0;side<2;side++){
       ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
       if(dir == dim-1){
-	ierr = IGABoundarySetUserIFunction(bnd,ZeroFluxResidual,&user);CHKERRQ(ierr);
-	ierr = IGABoundarySetUserIJacobian(bnd,ZeroFluxJacobian,&user);CHKERRQ(ierr);
-	if(side == 0){
-	  ierr = IGABoundarySetValue(bnd,0,user.Sin);CHKERRQ(ierr);
-	}
+        ierr = IGABoundarySetUserIFunction(bnd,ZeroFluxResidual,&user);CHKERRQ(ierr);
+        ierr = IGABoundarySetUserIJacobian(bnd,ZeroFluxJacobian,&user);CHKERRQ(ierr);
+        if(side == 0){
+          ierr = IGABoundarySetValue(bnd,0,user.Sin);CHKERRQ(ierr);
+        }
       }
     }
-  
+
   ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
   ierr = IGASetUserIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
 

File include/petiga.h

 PETSC_EXTERN PetscErrorCode IGAPointGetShapeFuns(IGAPoint point,PetscInt der,const PetscReal *shapefuns[]);
 
 PETSC_EXTERN PetscErrorCode IGAPointInterpolate(IGAPoint point,PetscInt ider,const PetscScalar U[],PetscScalar u[]);
-
 PETSC_EXTERN PetscErrorCode IGAPointFormPoint    (IGAPoint p,PetscReal x[]);
 PETSC_EXTERN PetscErrorCode IGAPointFormGradMap  (IGAPoint p,PetscReal map[],PetscReal inv[]);
 PETSC_EXTERN PetscErrorCode IGAPointFormShapeFuns(IGAPoint p,PetscInt der,PetscReal N[]);
 
 PETSC_EXTERN PetscErrorCode IGAPointGetWorkVec(IGAPoint point,PetscScalar *V[]);
 PETSC_EXTERN PetscErrorCode IGAPointGetWorkMat(IGAPoint point,PetscScalar *M[]);
-
 PETSC_EXTERN PetscErrorCode IGAPointAddArray(IGAPoint point,PetscInt n,const PetscScalar a[],PetscScalar A[]);
 PETSC_EXTERN PetscErrorCode IGAPointAddVec(IGAPoint point,const PetscScalar f[],PetscScalar F[]);
 PETSC_EXTERN PetscErrorCode IGAPointAddMat(IGAPoint point,const PetscScalar k[],PetscScalar K[]);