1. Lisandro Dalcin
  2. PetIGA

Commits

Lisandro Dalcin  committed 86361b0

Expose the hessian of the geometry map in IGAPoint

  • Participants
  • Parent commits 1144986
  • Branches default

Comments (0)

Files changed (9)

File include/petiga.h

View file
   PetscReal *detX;     /*   [nqp]                     */
   PetscReal *gradX[2]; /*0: [nqp][nsd][dim]           */
                        /*1: [nqp][dim][nsd]           */
+  PetscReal *hessX[2]; /*0: [nqp][nsd][dim][dim]      */
+                       /*1: [nqp][dim][nsd][nsd       */
   PetscReal *detS;     /*   [nqp]                     */
   PetscReal *normal;   /*   [nqp][dim]                */
 
   PetscReal *detX;     /*   [1] */
   PetscReal *gradX[2]; /*0: [nsd][dim] */
                        /*1: [dim][nsd] */
+  PetscReal *hessX[2]; /*0: [nsd][dim][dim] */
+                       /*1: [dim][nsd][nsd] */
   PetscReal *detS;     /*   [1] */
   PetscReal *normal;   /*   [dim] */
 

File src/petiga1d.f90

View file
      nqp,nen,X,                  &
      M0,M1,M2,M3,                &
      N0,N1,N2,N3,                &
-     DetF,F,G)                   &
+     dX,G0,G1,H0,H1)             &
   bind(C, name="IGA_ShapeFuns_1D")
   use PetIGA
   implicit none
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: DetF(nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: F(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
   integer(kind=IGA_INTEGER_KIND)  :: q
   do q=1,nqp
      call GeometryMap(&
           nen,X,&
           M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
           N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          DetF(q),F(:,:,q),G(:,:,q))
+          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
   end do
 contains
 include 'petigageo.f90.in'

File src/petiga2d.f90

View file
      nqp,nen,X,                  &
      M0,M1,M2,M3,                &
      N0,N1,N2,N3,                &
-     DetF,F,G)                   &
+     dX,G0,G1,H0,H1)             &
   bind(C, name="IGA_ShapeFuns_2D")
   use PetIGA
   implicit none
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: DetF(nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: F(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
   integer(kind=IGA_INTEGER_KIND)  :: q
   do q=1,nqp
      call GeometryMap(&
           nen,X,&
           M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
           N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          DetF(q),F(:,:,q),G(:,:,q))
+          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
   end do
 contains
 include 'petigageo.f90.in'

File src/petiga3d.f90

View file
      nqp,nen,X,                  &
      M0,M1,M2,M3,                &
      N0,N1,N2,N3,                &
-     DetF,F,G)                   &
+     dX,G0,G1,H0,H1)             &
   bind(C, name="IGA_ShapeFuns_3D")
   use PetIGA
   implicit none
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
   real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: DetF(nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: F(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
   integer(kind=IGA_INTEGER_KIND)  :: q
   do q=1,nqp
      call GeometryMap(&
           nen,X,&
           M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
           N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          DetF(q),F(:,:,q),G(:,:,q))
+          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
   end do
 contains
 include 'petigageo.f90.in'

File src/petigaelem.c

View file
   ierr = PetscFree(element->detX);CHKERRQ(ierr);
   ierr = PetscFree(element->gradX[0]);CHKERRQ(ierr);
   ierr = PetscFree(element->gradX[1]);CHKERRQ(ierr);
+  ierr = PetscFree(element->hessX[0]);CHKERRQ(ierr);
+  ierr = PetscFree(element->hessX[1]);CHKERRQ(ierr);
   ierr = PetscFree(element->detS);CHKERRQ(ierr);
   ierr = PetscFree(element->normal);CHKERRQ(ierr);
 
     ierr = PetscMalloc1(nqp,PetscReal,&element->detX);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*dim*dim,PetscReal,&element->gradX[0]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*dim*dim,PetscReal,&element->gradX[1]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*dim*dim*dim,PetscReal,&element->hessX[0]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nqp*dim*dim*dim,PetscReal,&element->hessX[1]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp,PetscReal,&element->detS);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*dim,PetscReal,&element->normal);CHKERRQ(ierr);
 
     ierr = PetscMemzero(element->detX,    sizeof(PetscReal)*nqp);CHKERRQ(ierr);
     ierr = PetscMemzero(element->gradX[0],sizeof(PetscReal)*nqp*dim*dim);CHKERRQ(ierr);
     ierr = PetscMemzero(element->gradX[1],sizeof(PetscReal)*nqp*dim*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(element->hessX[0],sizeof(PetscReal)*nqp*dim*dim*dim);CHKERRQ(ierr);
+    ierr = PetscMemzero(element->hessX[1],sizeof(PetscReal)*nqp*dim*dim*dim);CHKERRQ(ierr);
     ierr = PetscMemzero(element->detS,    sizeof(PetscReal)*nqp);CHKERRQ(ierr);
     ierr = PetscMemzero(element->normal,  sizeof(PetscReal)*nqp*dim);CHKERRQ(ierr);
     /* */
     ierr = PetscMemzero(element->shape[3],sizeof(PetscReal)*nqp*nen*dim*dim*dim);CHKERRQ(ierr);
     /* */
     for (q=0; q<nqp; q++) {
-      PetscReal *F = &element->gradX[0][q*dim*dim];
-      PetscReal *G = &element->gradX[1][q*dim*dim];
+      PetscReal *G0 = &element->gradX[0][q*dim*dim];
+      PetscReal *G1 = &element->gradX[1][q*dim*dim];
       element->detX[q] = 1.0;
       for (i=0; i<dim; i++)
-        F[i*(dim+1)] = G[i*(dim+1)] = 1.0;
+        G0[i*(dim+1)] = G1[i*(dim+1)] = 1.0;
     }
   }
   PetscFunctionReturn(0);
   point->detX     += 1;
   point->gradX[0] += dim*dim;
   point->gradX[1] += dim*dim;
+  point->hessX[0] += dim*dim*dim;
+  point->hessX[1] += dim*dim*dim;
   point->detS     += 1;
   point->normal   += dim;
 
   point->detX     = element->detX;
   point->gradX[0] = element->gradX[0];
   point->gradX[1] = element->gradX[1];
+  point->hessX[0] = element->hessX[0];
+  point->hessX[1] = element->hessX[1];
   point->detS     = element->detS;
   point->normal   = element->normal;
 
 extern void IGA_ShapeFuns_1D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                              const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                              PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-                             PetscReal[],PetscReal[],PetscReal[]);
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 extern void IGA_ShapeFuns_2D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                              const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                              PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-                             PetscReal[],PetscReal[],PetscReal[]);
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 extern void IGA_ShapeFuns_3D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                              const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                              PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-                             PetscReal[],PetscReal[],PetscReal[]);
+                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 EXTERN_C_END
 
 #define IGA_Quadrature_ARGS(ID,BD,i) \
     PetscReal **M = element->basis;
     PetscReal **N = element->shape;
     PetscReal *J  = element->detX;
-    PetscReal *F1 = element->gradX[0];
+    PetscReal *G0 = element->gradX[0];
     PetscReal *G1 = element->gradX[1];
+    PetscReal *H0 = element->hessX[0];
+    PetscReal *H1 = element->hessX[1];
     switch (element->dim) {
     case 3: IGA_ShapeFuns_3D(ord,nqp,nen,X,
                              M[0],M[1],M[2],M[3],
                              N[0],N[1],N[2],N[3],
-                             J,F1,G1); break;
+                             J,G0,G1,H0,H1); break;
     case 2: IGA_ShapeFuns_2D(ord,nqp,nen,X,
                              M[0],M[1],M[2],M[3],
                              N[0],N[1],N[2],N[3],
-                             J,F1,G1); break;
+                             J,G0,G1,H0,H1); break;
     case 1: IGA_ShapeFuns_1D(ord,nqp,nen,X,
                              M[0],M[1],M[2],M[3],
                              N[0],N[1],N[2],N[3],
-                             J,F1,G1); break;
+                             J,G0,G1,H0,H1); break;
     }
     for (q=0; q<nqp; q++)
       element->detJac[q] *= J[q];
     PetscReal **M = element->basis;
     PetscReal **N = element->shape;
     PetscReal *J  = element->detX;
-    PetscReal *F1 = element->gradX[0];
+    PetscReal *G0 = element->gradX[0];
     PetscReal *G1 = element->gradX[1];
+    PetscReal *H0 = element->hessX[0];
+    PetscReal *H1 = element->hessX[1];
     PetscReal *S  = element->detS;
     PetscReal *n  = element->normal;
     switch (dim) {
     case 3: IGA_ShapeFuns_3D(ord,nqp,nen,X,
                              M[0],M[1],M[2],M[3],
                              N[0],N[1],N[2],N[3],
-                             J,F1,G1); break;
+                             J,G0,G1,H0,H1); break;
     case 2: IGA_ShapeFuns_2D(ord,nqp,nen,X,
                              M[0],M[1],M[2],M[3],
                              N[0],N[1],N[2],N[3],
-                             J,F1,G1); break;
+                             J,G0,G1,H0,H1); break;
     case 1: IGA_ShapeFuns_1D(ord,nqp,nen,X,
                              M[0],M[1],M[2],M[3],
                              N[0],N[1],N[2],N[3],
-                             J,F1,G1); break;
+                             J,G0,G1,H0,H1); break;
     }
     for (q=0; q<nqp; q++)
-      IGA_GetNormal(dim,dir,side,&F1[q*dim*dim],&S[q],&n[q*dim]);
+      IGA_GetNormal(dim,dir,side,&G0[q*dim*dim],&S[q],&n[q*dim]);
     for (q=0; q<nqp; q++)
       element->detJac[q] *= S[q];
   }

File src/petigaftn.F90

View file
      type(C_PTR) :: basis(0:3)
      type(C_PTR) :: detX
      type(C_PTR) :: gradX(0:1)
+     type(C_PTR) :: hessX(0:1)
      type(C_PTR) :: detS
      type(C_PTR) :: normal
      type(C_PTR) :: shape(0:3)

File src/petigageo.f90.in

View file
      nen,X,&
      N0,N1,N2,N3,&
      R0,R1,R2,R3,&
-     DetG,Grad,InvG)
+     dX,X1,E1,X2,E2)
   use PetIGA
   implicit none
   !integer(kind=IGA_INTEGER_KIND),parameter   :: dim = 1,2,3
   real   (kind=IGA_REAL_KIND   ), intent(out) :: R1(        dim,nen)
   real   (kind=IGA_REAL_KIND   ), intent(out) :: R2(    dim,dim,nen)
   real   (kind=IGA_REAL_KIND   ), intent(out) :: R3(dim,dim,dim,nen)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: DetG
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: Grad(dim,dim)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: InvG(dim,dim)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: dX
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: X1(dim,dim)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: E1(dim,dim)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: X2(dim,dim,dim)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: E2(dim,dim,dim)
 
-  integer(kind=IGA_INTEGER_KIND)  :: idx
+  real   (kind=IGA_REAL_KIND   )  :: X3(dim,dim,dim,dim)
+  real   (kind=IGA_REAL_KIND   )  :: E3(dim,dim,dim,dim)
+  integer(kind=IGA_INTEGER_KIND)  :: node
   integer(kind=IGA_INTEGER_KIND)  :: i, j, k, l
   integer(kind=IGA_INTEGER_KIND)  :: a, b, c, d
-  real   (kind=IGA_REAL_KIND   )  :: X1(dim,dim)
-  real   (kind=IGA_REAL_KIND   )  :: X2(dim,dim,dim)
-  real   (kind=IGA_REAL_KIND   )  :: X3(dim,dim,dim,dim)
-  real   (kind=IGA_REAL_KIND   )  :: E1(dim,dim)
-  real   (kind=IGA_REAL_KIND   )  :: E2(dim,dim,dim)
-  real   (kind=IGA_REAL_KIND   )  :: E3(dim,dim,dim,dim)
 
   ! gradient of the mapping
-  Grad = matmul(N1,transpose(X))
-  DetG = Determinant(dim,Grad)
-  InvG = Inverse(dim,DetG,Grad)
+  X1 = matmul(N1,transpose(X))
+  dX = Determinant(dim,X1)
+  E1 = Inverse(dim,dX,X1)
 
   ! 0th derivatives
   R0 = N0
 
   ! 1st derivatives
-  X1 = transpose(Grad)
-  E1 = transpose(InvG)
-  R1 = 0
-  do idx = 1,nen
+  if (order < 1) return
+  R1 = 0 ! matmul(E1,N1)
+  do node = 1,nen
      do i = 1,dim
         do a = 1,dim
-           R1(i,idx) = N1(a,idx)*E1(a,i) +  R1(i,idx)
+           R1(i,node) = R1(i,node) + N1(a,node)*E1(i,a)
         end do
      end do
   end do
   ! 2nd derivatives
   if (order < 2) return
   X2 = 0
-  do b = 1,dim
-     do a = 1,dim
-        do i = 1,dim
-           do idx = 1,nen
-              X2(i,a,b) = X(i,idx)*N2(a,b,idx) + X2(i,a,b)
+  do node = 1,nen
+     do i = 1,dim
+        do a = 1,dim
+           do b = 1,dim
+              X2(b,a,i) = X2(b,a,i) + X(i,node)*N2(b,a,node)
            end do
         end do
      end do
   end do
   E2 = 0
-  do j = 1,dim
-     do i = 1,dim
-        do c = 1,dim
-           do b = 1,dim
-              do a = 1,dim
-                 do k = 1,dim
-                    E2(c,i,j) = - X2(k,a,b)*E1(a,i)*E1(b,j)*E1(c,k) + E2(c,i,j)
+  do i = 1,dim
+     do j = 1,dim
+        do k = 1,dim
+           do a = 1,dim
+              do b = 1,dim
+                 do c = 1,dim
+                    E2(j,i,c) = E2(j,i,c) - X2(b,a,k)*E1(i,a)*E1(j,b)*E1(k,c)
                  end do
               end do
            end do
      end do
   end do
   R2 = 0
-  do idx = 1,nen
-     do j = 1,dim
-        do i = 1,dim
-           do b = 1,dim
-              do a = 1,dim
-                 R2(i,j,idx) = N2(a,b,idx)*E1(a,i)*E1(b,j) + R2(i,j,idx)
+  do node = 1,nen
+     do i = 1,dim
+        do j = 1,dim
+           do a = 1,dim
+              do b = 1,dim
+                 R2(j,i,node) = R2(j,i,node) + N2(b,a,node)*E1(i,a)*E1(j,b)
               end do
-              R2(i,j,idx) = N1(b,idx)*E2(b,i,j) + R2(i,j,idx)
+              R2(j,i,node) = R2(j,i,node) + N1(a,node)*E2(j,i,a)
            end do
         end do
      end do
   ! 3rd derivatives
   if (order < 3) return
   X3 = 0
-  do c = 1,dim
-     do b = 1,dim
+  do node = 1,nen
+     do i = 1,dim
         do a = 1,dim
-           do i = 1,dim
-              do idx = 1,nen
-                 X3(i,a,b,c) = X(i,idx)*N3(a,b,c,idx) + X3(i,a,b,c)
+           do b = 1,dim
+              do c = 1,dim
+                 X3(c,b,a,i) = X3(c,b,a,i) &
+                      +X(i,node)*N3(c,b,a,node)
               end do
            end do
         end do
      end do
   end do
   E3 = 0
-  do k = 1,dim
-     do j = 1,dim
-        do i = 1,dim
-           do d = 1,dim
+  do d = 1,dim
+     do i = 1,dim
+        do j = 1,dim
+           do k = 1,dim
               do a = 1,dim
                  do b = 1,dim
                     do l = 1,dim
                        do c = 1,dim
-                          E3(d,i,j,k) = - X3(l,c,b,a)*E1(c,i)*E1(b,j)*E1(a,k)*E1(d,l) + E3(d,i,j,k)
+                          E3(k,j,i,d) = E3(k,j,i,d) &
+                               -X3(c,b,a,l)*E1(i,a)*E1(j,b)*E1(k,c)*E1(l,d)
                        end do
-                       E3(d,i,j,k) = - X2(l,b,a)*( &
-                            + E1(a,j)*E2(b,i,k) &
-                            + E1(a,k)*E2(b,i,j) &
-                            + E1(b,i)*E2(a,j,k) )*E1(d,l) + E3(d,i,j,k)
+                       E3(k,j,i,d) = E3(k,j,i,d) &
+                            -X2(b,a,l)*( &
+                            +E1(i,a)*E2(k,j,b) &
+                            +E1(j,b)*E2(k,i,a) &
+                            +E1(k,b)*E2(j,i,a))*E1(l,d)
                     end do
                  end do
               end do
      end do
   end do
   R3 = 0
-  do idx = 1,nen
-     do k = 1,dim
+  do node = 1,nen
+     do i = 1,dim
         do j = 1,dim
-           do i = 1,dim
+           do k = 1,dim
               do a = 1,dim
                  do b = 1,dim
                     do c = 1,dim
-                       R3(i,j,k,idx) = N3(c,b,a,idx)*E1(c,i)*E1(b,j)*E1(a,k) + R3(i,j,k,idx)
+                       R3(k,j,i,node) = R3(k,j,i,node) &
+                            +N3(c,b,a,node)*E1(i,a)*E1(j,b)*E1(k,c)
                     end do
-                    R3(i,j,k,idx) = N2(b,a,idx)*( &
-                         + E1(a,j)*E2(b,i,k) &
-                         + E1(a,k)*E2(b,i,j) &
-                         + E1(b,i)*E2(a,j,k) ) + R3(i,j,k,idx)
+                    R3(k,j,i,node) = R3(k,j,i,node) &
+                         +N2(b,a,node)*( &
+                         +E1(i,a)*E2(k,j,b) &
+                         +E1(j,b)*E2(k,i,a) &
+                         +E1(k,b)*E2(j,i,a))
                  end do
-                 R3(i,j,k,idx) = N1(a,idx)*E3(a,i,j,k) + R3(i,j,k,idx)
+                 R3(k,j,i,node) = R3(k,j,i,node) &
+                      +N1(a,node)*E3(k,j,i,a)
               end do
            end do
         end do

File src/petigapoint.c

View file
 extern void IGA_ShapeFuns_1D(PetscInt,PetscInt,PetscInt,const PetscReal[],
 			     const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
 			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[]);
+			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 extern void IGA_ShapeFuns_2D(PetscInt,PetscInt,PetscInt,const PetscReal[],
 			     const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
 			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[]);
+			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 extern void IGA_ShapeFuns_3D(PetscInt,PetscInt,PetscInt,const PetscReal[],
 			     const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
 			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],
-			     PetscReal[],PetscReal[],PetscReal[]);
+			     PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
 EXTERN_C_END
 
 #undef  __FUNCT__
       PetscReal **N = element->shape;
       PetscReal *dX = element->detX;
       PetscReal **gX = element->gradX;
+      PetscReal **hX = element->hessX;
       switch (element->dim) {
       case 3: IGA_ShapeFuns_3D(order,
 			       1,element->nen,
 			       element->geometryX,
 			       M[0],M[1],M[2],M[3],
 			       N[0],N[1],N[2],N[3],
-			       dX,gX[0],gX[1]); break;
+			       dX,gX[0],gX[1],hX[0],hX[1]); break;
       case 2: IGA_ShapeFuns_2D(order,
 			       1,element->nen,
 			       element->geometryX,
 			       M[0],M[1],M[2],M[3],
 			       N[0],N[1],N[2],N[3],
-			       dX,gX[0],gX[1]); break;
+			       dX,gX[0],gX[1],hX[0],hX[1]); break;
       case 1: IGA_ShapeFuns_1D(order,
 			       1,element->nen,
 			       element->geometryX,
 			       M[0],M[1],M[2],M[3],
 			       N[0],N[1],N[2],N[3],
-			       dX,gX[0],gX[1]); break;
+			       dX,gX[0],gX[1],hX[0],hX[1]); break;
       }
     }
 
     point->detX     = element->detX;
     point->gradX[0] = element->gradX[0];
     point->gradX[1] = element->gradX[1];
+    point->hessX[0] = element->hessX[0];
+    point->hessX[1] = element->hessX[1];
     point->shape[0] = element->shape[0];
     point->shape[1] = element->shape[1];
     point->shape[2] = element->shape[2];

File test/GeometryMap.c

View file
     AssertEQUAL(F[2][2], 2.0);
   }
   if (dim==2) {
-    PetscInt i,j,k;
-    PetscInt nen = p->nen;
+    PetscInt  i,j,k;
     PetscReal G[2][2];
     PetscReal H[2][2][2];
-    ierr = PetscMemcpy(&G[0][0],p->gradX[0],sizeof(G));CHKERRQ(ierr);
-    {
-      PetscInt  a;
-      PetscReal (*C)[2]    = (PetscReal(*)[2])    p->geometry;
-      PetscReal (*N)[2][2] = (PetscReal(*)[2][2]) p->basis[2];
-      ierr = PetscMemzero(&H[0][0][0],sizeof(H));CHKERRQ(ierr);
-      for (a=0;a<nen;a++)
-        for (k=0;k<dim;k++)
-          for (i=0;i<dim;i++)
-            for (j=0;j<dim;j++)
-              H[k][i][j] += C[a][k]*N[a][i][j];
-    }
+    ierr = PetscMemcpy(&G[0][0],   p->gradX[0],sizeof(G));CHKERRQ(ierr);
+    ierr = PetscMemcpy(&H[0][0][0],p->hessX[0],sizeof(H));CHKERRQ(ierr);
     for (k=0;k<dim;k++)
       for (i=0;i<dim;i++)
         for (j=0;j<dim;j++)
   }
   if (dim==3) {
     PetscInt i,j,k;
-    PetscInt nen = p->nen;
     PetscReal G[3][3];
     PetscReal H[3][3][3];
-    ierr = PetscMemcpy(&G[0][0],p->gradX[0],sizeof(G));CHKERRQ(ierr);
-    {
-      PetscInt  a;
-      PetscReal (*C)[3]    = (PetscReal(*)[3])    p->geometry;
-      PetscReal (*N)[3][3] = (PetscReal(*)[3][3]) p->basis[2];
-      ierr = PetscMemzero(&H[0][0][0],sizeof(H));CHKERRQ(ierr);
-      for (a=0;a<nen;a++)
-        for (k=0;k<dim;k++)
-          for (i=0;i<dim;i++)
-            for (j=0;j<dim;j++)
-              H[k][i][j] += C[a][k]*N[a][i][j];
-    }
+    ierr = PetscMemcpy(&G[0][0],   p->gradX[0],sizeof(G));CHKERRQ(ierr);
+    ierr = PetscMemcpy(&H[0][0][0],p->hessX[0],sizeof(H));CHKERRQ(ierr);
     for (k=0;k<dim;k++)
       for (i=0;i<dim;i++)
         for (j=0;j<dim;j++)
       AssertEQUAL(kappa, 1/radius);
     }
   }
+  if (dim==2) {
+    PetscInt  nen = p->nen;
+    PetscInt  a,k,i,j;
+    PetscReal (*C)[2] = (PetscReal(*)[2])    p->geometry;
+    PetscReal (*N1)[2] = (PetscReal(*)[2]) p->shape[1];
+    PetscReal (*N2)[2][2] = (PetscReal(*)[2][2]) p->shape[2];
+    PetscReal G[2][2],H[2][2][2];
+    ierr = PetscMemzero(&G[0][0],sizeof(G));CHKERRQ(ierr);
+    ierr = PetscMemzero(&H[0][0][0],sizeof(H));CHKERRQ(ierr);
+    for (a=0;a<nen;a++)
+      for (k=0;k<dim;k++)
+        for (i=0;i<dim;i++)
+          G[k][i] += C[a][k]*N1[a][i];
+    for (a=0;a<nen;a++)
+      for (k=0;k<dim;k++)
+        for (i=0;i<dim;i++)
+          for (j=0;j<dim;j++)
+            H[k][i][j] += C[a][k]*N2[a][i][j];
+    for (i=0;i<dim;i++)
+      for (j=0;j<dim;j++)
+        if (i==j)
+          AssertEQUAL(G[i][j], 1.0);
+        else
+          AssertEQUAL(G[i][j], 0.0);
+    for (i=0;i<dim;i++)
+      for (j=0;j<dim;j++)
+        for (k=0;k<dim;k++)
+          AssertEQUAL(H[i][j][k], 0.0);
+  }
+  if (dim==3) {
+    PetscInt  nen = p->nen;
+    PetscInt  a,k,i,j;
+    PetscReal (*C)[3] = (PetscReal(*)[3])    p->geometry;
+    PetscReal (*N1)[3] = (PetscReal(*)[3]) p->shape[1];
+    PetscReal (*N2)[3][3] = (PetscReal(*)[3][3]) p->shape[2];
+    PetscReal G[3][3],H[3][3][3];
+    ierr = PetscMemzero(&G[0][0],sizeof(G));CHKERRQ(ierr);
+    ierr = PetscMemzero(&H[0][0][0],sizeof(H));CHKERRQ(ierr);
+    for (a=0;a<nen;a++)
+      for (k=0;k<dim;k++)
+        for (i=0;i<dim;i++)
+          G[k][i] += C[a][k]*N1[a][i];
+    for (k=0;k<dim;k++)
+      for (i=0;i<dim;i++)
+        if (k==i) AssertEQUAL(G[k][i], 1.0);
+        else      AssertEQUAL(G[k][i], 0.0);
+    for (a=0;a<nen;a++)
+      for (k=0;k<dim;k++)
+        for (i=0;i<dim;i++)
+          for (j=0;j<dim;j++)
+            H[k][i][j] += C[a][k]*N2[a][i][j];
+    for (k=0;k<dim;k++)
+      for (i=0;i<dim;i++)
+        for (j=0;j<dim;j++)
+          AssertEQUAL(H[k][i][j], 0.0);
+  }
   PetscFunctionReturn(0);
 }