Lisandro Dalcin avatar Lisandro Dalcin committed 340aa31

Fix IGAPointFormInvGradGeomMap() for dim!=nsd

Comments (0)

Files changed (3)

 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 IGAPointFormGeomMap(IGAPoint p,PetscReal x[]);
+PETSC_EXTERN PetscErrorCode IGAPointFormGradGeomMap(IGAPoint p,PetscReal F[]);
+PETSC_EXTERN PetscErrorCode IGAPointFormInvGradGeomMap(IGAPoint p,PetscReal G[]);
+
 PETSC_EXTERN PetscErrorCode IGAPointFormValue(IGAPoint p,const PetscScalar U[],PetscScalar u[]);
 PETSC_EXTERN PetscErrorCode IGAPointFormGrad (IGAPoint p,const PetscScalar U[],PetscScalar u[]);
 PETSC_EXTERN PetscErrorCode IGAPointFormHess (IGAPoint p,const PetscScalar U[],PetscScalar u[]);

src/petigapoint.c

     }
     for (a=0; a<dim; a++)
       for (i=0; i<nsd; i++)
-        G[a*dim+i] /= L[a];
+        G[a*nsd+i] /= L[a];
   } else {
     PetscInt i,dim = p->dim;
     const PetscReal *L = p->scale;

src/petigaval.F90

   real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
   real   (kind=IGA_REAL_KIND   ), intent(in)       :: C(nsd,nen)
   real   (kind=IGA_REAL_KIND   ), intent(out)      :: G(nsd,dim)
-  real   (kind=IGA_REAL_KIND   )  :: F(dim,nsd)
-  real   (kind=IGA_REAL_KIND   )  :: M(nsd,nsd), invM(nsd,nsd)
-  F = matmul(N,transpose(C))
+  real   (kind=IGA_REAL_KIND   )  :: F(nsd,dim)
+  real   (kind=IGA_REAL_KIND   )  :: M(dim,dim), detM, invM(dim,dim)
+  F = matmul(C,transpose(N))
   M = matmul(transpose(F),F)
-  invM = Inverse(nsd,Determinant(nsd,M),M)
-  G = matmul(invM,transpose(F))
+  detM = Determinant(dim,M)
+  invM = Inverse(dim,detM,M)
+  G = transpose(matmul(invM,transpose(F)))
 contains
 include 'petigainv.f90.in'
 end subroutine IGA_GetInvGradGeomMap
   bind(C, name="IGA_GetValue")
   use PetIGA
   implicit none
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof9*
   real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(nen)
   scalar (kind=IGA_SCALAR_KIND ), intent(in)       :: U(dof,nen)
   scalar (kind=IGA_SCALAR_KIND ), intent(out)      :: V(dof)
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.