# PetIGA

committed 5e19e94

Rework geometric mapping on manifolds (dim<nsd)

# src/makefile

` `
` SOURCEH  = ../include/petiga.h petigapc.h petigagrid.h petigapart.h`
` SOURCEC  = petiga.c petigareg.c petigaftnc.c petigaaxis.c petigarule.c petigabasis.c petigabound.c petigaelem.c petigapoint.c petigavec.c petigamat.c petigascl.c petigapcb.c petigapce.c petigaksp.c petigasnes.c petigats.c petigaio.c petigagrid.c petigapart.c petigacol.c`
`-SOURCEF1 = petigaftn.F90 petigaint.F90`
`+SOURCEF1 = petigaftn.F90 petigaval.F90`
` SOURCEF2 = petigabsp.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90`
` SOURCEF  = \${SOURCEF1} \${SOURCEF2}`
` OBJSC    = \${SOURCEC:.c=.o}`

# src/petigageo.f90.in

` `
` end subroutine GeometryMap`
` `
`-pure function Determinant(dim,A) result (detA)`
`-  implicit none`
`-  integer(kind=IGA_INT ), intent(in) :: dim`
`-  real   (kind=IGA_REAL), intent(in) :: A(dim,dim)`
`-  real   (kind=IGA_REAL) :: detA`
`-  select case (dim)`
`-  case (1)`
`-     detA = A(1,1)`
`-  case (2)`
`-     detA = + A(1,1)*A(2,2) - A(2,1)*A(1,2)`
`-  case (3)`
`-     detA = + A(1,1) * ( A(2,2)*A(3,3) - A(3,2)*A(2,3) ) &`
`-            - A(2,1) * ( A(1,2)*A(3,3) - A(3,2)*A(1,3) ) &`
`-            + A(3,1) * ( A(1,2)*A(2,3) - A(2,2)*A(1,3) )`
`-  case default`
`-     detA = 0`
`-  end select`
`-end function Determinant`
`-`
`-pure function Inverse(dim,detA,A) result (invA)`
`-  implicit none`
`-  integer(kind=IGA_INT ), intent(in) :: dim`
`-  real   (kind=IGA_REAL), intent(in) :: detA`
`-  real   (kind=IGA_REAL), intent(in) :: A(dim,dim)`
`-  real   (kind=IGA_REAL)             :: invA(dim,dim)`
`-  select case (dim)`
`-  case (1)`
`-     invA = 1.0/detA`
`-  case (2)`
`-     invA(1,1) = + A(2,2)`
`-     invA(2,1) = - A(2,1)`
`-     invA(1,2) = - A(1,2)`
`-     invA(2,2) = + A(1,1)`
`-     invA = invA/detA`
`-  case (3)`
`-     invA(1,1) = + A(2,2)*A(3,3) - A(2,3)*A(3,2)`
`-     invA(2,1) = - A(2,1)*A(3,3) + A(2,3)*A(3,1)`
`-     invA(3,1) = + A(2,1)*A(3,2) - A(2,2)*A(3,1)`
`-     invA(1,2) = - A(1,2)*A(3,3) + A(1,3)*A(3,2)`
`-     invA(2,2) = + A(1,1)*A(3,3) - A(1,3)*A(3,1)`
`-     invA(3,2) = - A(1,1)*A(3,2) + A(1,2)*A(3,1)`
`-     invA(1,3) = + A(1,2)*A(2,3) - A(1,3)*A(2,2)`
`-     invA(2,3) = - A(1,1)*A(2,3) + A(1,3)*A(2,1)`
`-     invA(3,3) = + A(1,1)*A(2,2) - A(1,2)*A(2,1)`
`-     invA = invA/detA`
`-  case default`
`-     invA = 0`
`-  end select`
`-end function Inverse`
`+include 'petigainv.f90.in'`

# src/petigaint.F90

`-! -*- f90 -*-`
`-`
`-#include "petscconf.h"`
`-#if defined(PETSC_USE_COMPLEX)`
`-#define scalar complex`
`-#else`
`-#define scalar real`
`-#endif`
`-`
`-subroutine IGA_GetPoint(nen,dim,N,C,X) &`
`-  bind(C, name="IGA_GetPoint")`
`-  use PetIGA`
`-  implicit none`
`-  integer(kind=IGA_INT ), intent(in),value :: nen,dim`
`-  real   (kind=IGA_REAL), intent(in)       :: N(nen)`
`-  real   (kind=IGA_REAL), intent(in)       :: C(dim,nen)`
`-  real   (kind=IGA_REAL), intent(out)      :: X(dim)`
`-  integer(kind=IGA_INT )  :: a`
`-  ! X = matmul(C,N)`
`-  X = 0`
`-  do a = 1, nen`
`-     X = X + N(a) * C(:,a)`
`-  end do`
`-end subroutine IGA_GetPoint`
`-`
`-subroutine IGA_GetValue(nen,dof,N,U,V) &`
`-  bind(C, name="IGA_GetValue")`
`-  use PetIGA`
`-  implicit none`
`-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof`
`-  real   (kind=IGA_REAL  ), intent(in)       :: N(nen)`
`-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dof)`
`-  integer(kind=IGA_INT   )  :: a, i`
`-  ! V = matmul(N,transpose(U))`
`-  V = 0`
`-  do a = 1, nen`
`-     V = V + N(a) * U(:,a)`
`-  end do`
`-end subroutine IGA_GetValue`
`-`
`-subroutine IGA_GetGrad(nen,dof,dim,N,U,V) &`
`-  bind(C, name="IGA_GetGrad")`
`-  use PetIGA`
`-  implicit none`
`-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim`
`-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim,nen)`
`-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim,dof)`
`-  integer(kind=IGA_INT   )  :: a, c`
`-  ! V = matmul(N,transpose(U))`
`-  V = 0`
`-  do a = 1, nen`
`-     do c = 1, dof`
`-        V(:,c) = V(:,c) + N(:,a) * U(c,a)`
`-     end do`
`-  end do`
`-end subroutine IGA_GetGrad`
`-`
`-subroutine IGA_GetHess(nen,dof,dim,N,U,V) &`
`-  bind(C, name="IGA_GetHess")`
`-  use PetIGA`
`-  implicit none`
`-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim`
`-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim*dim,nen)`
`-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim*dim,dof)`
`-  integer(kind=IGA_INT   )  :: a, i`
`-  ! V = matmul(N,transpose(U))`
`-  V = 0`
`-  do a = 1, nen`
`-     do i = 1, dof`
`-        V(:,i) = V(:,i) + N(:,a) * U(i,a)`
`-     end do`
`-  end do`
`-end subroutine IGA_GetHess`
`-`
`-subroutine IGA_GetDel2(nen,dof,dim,N,U,V) &`
`-  bind(C, name="IGA_GetDel2")`
`-  use PetIGA`
`-  implicit none`
`-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim`
`-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim,dim,nen)`
`-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dof)`
`-  integer(kind=IGA_INT   )  :: a, c, i`
`-  V = 0`
`-  do a = 1, nen`
`-     do c = 1, dof`
`-        do i = 1, dim`
`-           V(c) = V(c) + N(i,i,a) * U(c,a)`
`-        end do`
`-     end do`
`-  end do`
`-end subroutine IGA_GetDel2`
`-`
`-subroutine IGA_GetDer3(nen,dof,dim,N,U,V) &`
`-     bind(C, name="IGA_GetDer3")`
`-  use PetIGA`
`-  implicit none`
`-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim`
`-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim*dim*dim,nen)`
`-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim*dim*dim,dof)`
`-  integer(kind=IGA_INT   )  :: a, i`
`-  ! V = matmul(N,transpose(U))`
`-  V = 0`
`-  do a = 1, nen`
`-     do i = 1, dof`
`-        V(:,i) = V(:,i) + N(:,a) * U(i,a)`
`-     end do`
`-  end do`
`-end subroutine IGA_GetDer3`
`-`
`-!subroutine IGA_GetDerivative(nen,dof,dim,der,N,U,V) &`
`-!  bind(C, name="IGA_GetDerivative")`
`-!  use PetIGA`
`-!  implicit none`
`-!  integer(kind=IGA_INT   ), intent(in),value :: nen,dof`
`-!  integer(kind=IGA_INT   ), intent(in),value :: dim,der`
`-!  real   (kind=IGA_REAL  ), intent(in)       :: N(dim**der,nen)`
`-!  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`-!  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim**der,dof)`
`-!  integer(kind=IGA_INT   )  :: a, i`
`-!  ! V = matmul(N,transpose(U))`
`-!  V = 0`
`-!  do a = 1, nen`
`-!     do i = 1, dof`
`-!        V(:,i) = V(:,i) + N(:,a) * U(i,a)`
`-!     end do`
`-!  end do`
`-!end subroutine IGA_GetDerivative`
`-`
`-subroutine IGA_Interpolate(nen,dof,dim,der,N,U,V) &`
`-  bind(C, name="IGA_Interpolate")`
`-  use PetIGA`
`-  implicit none`
`-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof`
`-  integer(kind=IGA_INT   ), intent(in),value :: dim,der`
`-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim**der,nen)`
`-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim**der,dof)`
`-  integer(kind=IGA_INT   )  :: a, i`
`-  ! V = matmul(N,transpose(U))`
`-  V = 0`
`-  do a = 1, nen`
`-     do i = 1, dof`
`-        V(:,i) = V(:,i) + N(:,a) * U(i,a)`
`-     end do`
`-  end do`
`-end subroutine IGA_Interpolate`

# src/petigainv.f90.in

`+! -*- f90 -*-`
`+`
`+pure function Determinant(dim,A) result (detA)`
`+  use PetIGA; implicit none`
`+  integer(kind=IGA_INT ), intent(in) :: dim`
`+  real   (kind=IGA_REAL), intent(in) :: A(dim,dim)`
`+  real   (kind=IGA_REAL)  :: detA`
`+  select case (dim)`
`+  case (1)`
`+     detA = A(1,1)`
`+  case (2)`
`+     detA = + A(1,1)*A(2,2) - A(2,1)*A(1,2)`
`+  case (3)`
`+     detA = + A(1,1) * ( A(2,2)*A(3,3) - A(3,2)*A(2,3) ) &`
`+            - A(2,1) * ( A(1,2)*A(3,3) - A(3,2)*A(1,3) ) &`
`+            + A(3,1) * ( A(1,2)*A(2,3) - A(2,2)*A(1,3) )`
`+  case default`
`+     detA = 0`
`+  end select`
`+end function Determinant`
`+`
`+pure function Inverse(dim,detA,A) result (invA)`
`+  use PetIGA; implicit none`
`+  integer(kind=IGA_INT ), intent(in) :: dim`
`+  real   (kind=IGA_REAL), intent(in) :: detA`
`+  real   (kind=IGA_REAL), intent(in) :: A(dim,dim)`
`+  real   (kind=IGA_REAL)  :: invA(dim,dim)`
`+  select case (dim)`
`+  case (1)`
`+     invA = 1.0/detA`
`+  case (2)`
`+     invA(1,1) = + A(2,2)`
`+     invA(2,1) = - A(2,1)`
`+     invA(1,2) = - A(1,2)`
`+     invA(2,2) = + A(1,1)`
`+     invA = invA/detA`
`+  case (3)`
`+     invA(1,1) = + A(2,2)*A(3,3) - A(2,3)*A(3,2)`
`+     invA(2,1) = - A(2,1)*A(3,3) + A(2,3)*A(3,1)`
`+     invA(3,1) = + A(2,1)*A(3,2) - A(2,2)*A(3,1)`
`+     invA(1,2) = - A(1,2)*A(3,3) + A(1,3)*A(3,2)`
`+     invA(2,2) = + A(1,1)*A(3,3) - A(1,3)*A(3,1)`
`+     invA(3,2) = - A(1,1)*A(3,2) + A(1,2)*A(3,1)`
`+     invA(1,3) = + A(1,2)*A(2,3) - A(1,3)*A(2,2)`
`+     invA(2,3) = - A(1,1)*A(2,3) + A(1,3)*A(2,1)`
`+     invA(3,3) = + A(1,1)*A(2,2) - A(1,2)*A(2,1)`
`+     invA = invA/detA`
`+  case default`
`+     invA = 0`
`+  end select`
`+end function Inverse`

# src/petigapoint.c

` }`
` `
` EXTERN_C_BEGIN`
`-extern void IGA_GetPoint(PetscInt nen,PetscInt dim,const PetscReal N[],`
`-                         const PetscReal Xw[],PetscReal x[]);`
`+extern void IGA_GetGeomMap (PetscInt nen,PetscInt nsd,`
`+                            const PetscReal N[],const PetscReal C[],PetscReal X[]);`
`+extern void IGA_GetGradMap (PetscInt nen,PetscInt nsd,PetscInt dim,`
`+                            const PetscReal N[],const PetscReal C[],PetscReal F[]);`
`+extern void IGA_GetGradMapI(PetscInt nen,PetscInt nsd,PetscInt dim,`
`+                            const PetscReal N[],const PetscReal C[],PetscReal G[]);`
`+EXTERN_C_END`
`+`
`+EXTERN_C_BEGIN`
` extern void IGA_GetValue(PetscInt nen,PetscInt dof,const PetscReal N[],`
`                          const PetscScalar U[],PetscScalar u[]);`
` extern void IGA_GetGrad (PetscInt nen,PetscInt dof,PetscInt dim,const PetscReal N[],`
`   PetscValidRealPointer(x,2);`
`   if (p->parent->geometry) {`
`     const PetscReal *X = p->parent->geometryX;`
`-    IGA_GetPoint(p->nen,p->dim,p->shape[0],X,x);`
`+    IGA_GetGeomMap(p->nen,p->nsd,p->shape[0],X,x);`
`   } else {`
`     PetscInt i,dim = p->dim;`
`     const PetscReal *X = p->point;`
`   if(F) PetscValidRealPointer(F,2);`
`   if(G) PetscValidRealPointer(G,3);`
`   if (p->parent->geometry) {`
`-    PetscInt i,a,dim = p->dim;`
`+    PetscInt a,dim = p->dim;`
`+    PetscInt i,nsd = p->nsd;`
`     const PetscReal *L = p->scale;`
`-    const PetscReal *X = p->gradX[0];`
`-    const PetscReal *Y = p->gradX[1];`
`-    if (F) {`
`-      for (i=0; i<dim; i++)`
`+    if (dim == nsd) {`
`+      if (F) {(void)PetscMemcpy(F,p->gradX[0],nsd*dim*sizeof(PetscReal));}`
`+      if (G) {(void)PetscMemcpy(G,p->gradX[1],dim*nsd*sizeof(PetscReal));}`
`+    } else {`
`+      const PetscReal *X = p->parent->geometryX;`
`+      if (F) {IGA_GetGradMap (p->nen,nsd,dim,p->basis[1],X,F);}`
`+      if (G) {IGA_GetGradMapI(p->nen,nsd,dim,p->basis[1],X,G);}`
`+    }`
`+    if (F)`
`+      for (i=0; i<nsd; i++)`
`         for (a=0; a<dim; a++)`
`-          F[i*dim+a] = X[i*dim+a]*L[a];`
`-    }`
`-    if (G) {`
`+          F[i*dim+a] *= L[a];`
`+    if (G)`
`       for (a=0; a<dim; a++)`
`-        for (i=0; i<dim; i++)`
`-          G[a*dim+i] = Y[a*dim+i]/L[a];`
`-    }`
`+        for (i=0; i<nsd; i++)`
`+          G[a*dim+i] /= L[a];`
`   } else {`
`     PetscInt i,dim = p->dim;`
`     const PetscReal *L = p->scale;`

# src/petigaval.F90

`+! -*- f90 -*-`
`+`
`+#include "petscconf.h"`
`+#if defined(PETSC_USE_COMPLEX)`
`+#define scalar complex`
`+#else`
`+#define scalar real`
`+#endif`
`+`
`+subroutine IGA_GetGeomMap(nen,nsd,N,C,X) &`
`+  bind(C, name="IGA_GetGeomMap")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT ), intent(in),value :: nen,nsd`
`+  real   (kind=IGA_REAL), intent(in)       :: N(nen)`
`+  real   (kind=IGA_REAL), intent(in)       :: C(nsd,nen)`
`+  real   (kind=IGA_REAL), intent(out)      :: X(nsd)`
`+  !integer(kind=IGA_INT )  :: a`
`+  !X = 0`
`+  !do a = 1, nen`
`+  !   X = X + N(a) * C(:,a)`
`+  !end do`
`+  X = matmul(C,N)`
`+end subroutine IGA_GetGeomMap`
`+`
`+subroutine IGA_GetGradMap(nen,nsd,dim,N,C,F) &`
`+  bind(C, name="IGA_GetGradMap")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT ), intent(in),value :: nen,nsd,dim`
`+  real   (kind=IGA_REAL), intent(in)       :: N(dim,nen)`
`+  real   (kind=IGA_REAL), intent(in)       :: C(nsd,nen)`
`+  real   (kind=IGA_REAL), intent(out)      :: F(dim,nsd)`
`+  F = matmul(N,transpose(C))`
`+end subroutine IGA_GetGradMap`
`+`
`+subroutine IGA_GetGradMapI(nen,nsd,dim,N,C,G) &`
`+  bind(C, name="IGA_GetGradMapI")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT ), intent(in),value :: nen,nsd,dim`
`+  real   (kind=IGA_REAL), intent(in)       :: N(dim,nen)`
`+  real   (kind=IGA_REAL), intent(in)       :: C(nsd,nen)`
`+  real   (kind=IGA_REAL), intent(out)      :: G(nsd,dim)`
`+  real   (kind=IGA_REAL)  :: F(dim,nsd)`
`+  real   (kind=IGA_REAL)  :: M(nsd,nsd), invM(nsd,nsd)`
`+  F = matmul(N,transpose(C))`
`+  M = matmul(transpose(F),F)`
`+  invM = Inverse(nsd,Determinant(nsd,M),M)`
`+  G = matmul(invM,transpose(F))`
`+contains`
`+include 'petigainv.f90.in'`
`+end subroutine IGA_GetGradMapI`
`+`
`+`
`+subroutine IGA_GetValue(nen,dof,N,U,V) &`
`+  bind(C, name="IGA_GetValue")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT   ), intent(in),value :: nen,dof`
`+  real   (kind=IGA_REAL  ), intent(in)       :: N(nen)`
`+  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`+  scalar (kind=IGA_SCALAR), intent(out)      :: V(dof)`
`+  integer(kind=IGA_INT   )  :: a, i`
`+  ! V = matmul(N,transpose(U))`
`+  V = 0`
`+  do a = 1, nen`
`+     V = V + N(a) * U(:,a)`
`+  end do`
`+end subroutine IGA_GetValue`
`+`
`+subroutine IGA_GetGrad(nen,dof,dim,N,U,V) &`
`+  bind(C, name="IGA_GetGrad")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim`
`+  real   (kind=IGA_REAL  ), intent(in)       :: N(dim,nen)`
`+  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`+  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim,dof)`
`+  integer(kind=IGA_INT   )  :: a, c`
`+  ! V = matmul(N,transpose(U))`
`+  V = 0`
`+  do a = 1, nen`
`+     do c = 1, dof`
`+        V(:,c) = V(:,c) + N(:,a) * U(c,a)`
`+     end do`
`+  end do`
`+end subroutine IGA_GetGrad`
`+`
`+subroutine IGA_GetHess(nen,dof,dim,N,U,V) &`
`+  bind(C, name="IGA_GetHess")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim`
`+  real   (kind=IGA_REAL  ), intent(in)       :: N(dim*dim,nen)`
`+  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`+  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim*dim,dof)`
`+  integer(kind=IGA_INT   )  :: a, i`
`+  ! V = matmul(N,transpose(U))`
`+  V = 0`
`+  do a = 1, nen`
`+     do i = 1, dof`
`+        V(:,i) = V(:,i) + N(:,a) * U(i,a)`
`+     end do`
`+  end do`
`+end subroutine IGA_GetHess`
`+`
`+subroutine IGA_GetDel2(nen,dof,dim,N,U,V) &`
`+  bind(C, name="IGA_GetDel2")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim`
`+  real   (kind=IGA_REAL  ), intent(in)       :: N(dim,dim,nen)`
`+  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`+  scalar (kind=IGA_SCALAR), intent(out)      :: V(dof)`
`+  integer(kind=IGA_INT   )  :: a, c, i`
`+  V = 0`
`+  do a = 1, nen`
`+     do c = 1, dof`
`+        do i = 1, dim`
`+           V(c) = V(c) + N(i,i,a) * U(c,a)`
`+        end do`
`+     end do`
`+  end do`
`+end subroutine IGA_GetDel2`
`+`
`+subroutine IGA_GetDer3(nen,dof,dim,N,U,V) &`
`+     bind(C, name="IGA_GetDer3")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim`
`+  real   (kind=IGA_REAL  ), intent(in)       :: N(dim*dim*dim,nen)`
`+  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`+  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim*dim*dim,dof)`
`+  integer(kind=IGA_INT   )  :: a, i`
`+  ! V = matmul(N,transpose(U))`
`+  V = 0`
`+  do a = 1, nen`
`+     do i = 1, dof`
`+        V(:,i) = V(:,i) + N(:,a) * U(i,a)`
`+     end do`
`+  end do`
`+end subroutine IGA_GetDer3`
`+`
`+!subroutine IGA_GetDerivative(nen,dof,dim,der,N,U,V) &`
`+!  bind(C, name="IGA_GetDerivative")`
`+!  use PetIGA`
`+!  implicit none`
`+!  integer(kind=IGA_INT   ), intent(in),value :: nen,dof`
`+!  integer(kind=IGA_INT   ), intent(in),value :: dim,der`
`+!  real   (kind=IGA_REAL  ), intent(in)       :: N(dim**der,nen)`
`+!  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`+!  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim**der,dof)`
`+!  integer(kind=IGA_INT   )  :: a, i`
`+!  ! V = matmul(N,transpose(U))`
`+!  V = 0`
`+!  do a = 1, nen`
`+!     do i = 1, dof`
`+!        V(:,i) = V(:,i) + N(:,a) * U(i,a)`
`+!     end do`
`+!  end do`
`+!end subroutine IGA_GetDerivative`
`+`
`+subroutine IGA_Interpolate(nen,dof,dim,der,N,U,V) &`
`+  bind(C, name="IGA_Interpolate")`
`+  use PetIGA`
`+  implicit none`
`+  integer(kind=IGA_INT   ), intent(in),value :: nen,dof`
`+  integer(kind=IGA_INT   ), intent(in),value :: dim,der`
`+  real   (kind=IGA_REAL  ), intent(in)       :: N(dim**der,nen)`
`+  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)`
`+  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim**der,dof)`
`+  integer(kind=IGA_INT   )  :: a, i`
`+  ! V = matmul(N,transpose(U))`
`+  V = 0`
`+  do a = 1, nen`
`+     do i = 1, dof`
`+        V(:,i) = V(:,i) + N(:,a) * U(i,a)`
`+     end do`
`+  end do`
`+end subroutine IGA_Interpolate`
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.