Commits

Lisandro Dalcin committed 4a4db86

A few changes to the handling of geometric mappings

Comments (0)

Files changed (8)

                        /*2: [nqp][nen][dim][dim]      */
                        /*3: [nqp][nen][dim][dim][dim] */
 
+  PetscReal *detX;     /*   [nqp]                     */
   PetscReal *gradX[2]; /*0: [nqp][nsd][dim]           */
                        /*1: [nqp][dim][nsd]           */
   PetscReal *shape[4]; /*0: [nqp][nen]                */
                        /*1: [nen][dim] */
                        /*2: [nen][dim][dim] */
                        /*3: [nen][dim][dim][dim] */
+
+  PetscReal *detX;     /*   [nqp] */
   PetscReal *gradX[2]; /*0: [nsd][dim] */
                        /*1: [dim][nsd] */
   PetscReal *shape[4]; /*0: [nen]  */
                         /*2: [nen][dim][dim]      */
                         /*3: [nen][dim][dim][dim] */
 
+  PetscReal detX;       /*   [1]                  */  
   PetscReal *gradX[2];  /*0: [nsd][dim]           */
                         /*1: [dim][nsd]           */
   PetscReal *shape[4];  /*0: [nen]                */
      nqp,nen,X,             &
      M0,M1,M2,M3,           &
      N0,N1,N2,N3,           &
-     DetJac,F,G)            &
+     DetF,F,G)              &
   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(inout) :: DetJac(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)
   integer(kind=IGA_INTEGER_KIND)  :: q
-  real   (kind=IGA_REAL_KIND)     :: DetF
   do q=1,nqp
      call GeometryMap(&
           order,&
           nen,X,&
           M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
           N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          DetF,F(:,:,q),G(:,:,q))
-     DetJac(q) = DetJac(q) * DetF
+          DetF(q),F(:,:,q),G(:,:,q))
   end do
 contains
 include 'petigageo.f90.in'
      nqp,nen,X,             &
      M0,M1,M2,M3,           &
      N0,N1,N2,N3,           &
-     DetJac,F,G)            &
+     DetF,F,G)              &
   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(inout) :: DetJac(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)
   integer(kind=IGA_INTEGER_KIND)  :: q
-  real   (kind=IGA_REAL_KIND)     :: DetF
   do q=1,nqp
      call GeometryMap(&
           order,&
           nen,X,&
           M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
           N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          DetF,F(:,:,q),G(:,:,q))
-     DetJac(q) = DetJac(q) * DetF
+          DetF(q),F(:,:,q),G(:,:,q))
   end do
 contains
 include 'petigageo.f90.in'
      nqp,nen,X,             &
      M0,M1,M2,M3,           &
      N0,N1,N2,N3,           &
-     DetJac,F,G)            &
+     DetF,F,G)              &
   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(inout) :: DetJac(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)
   integer(kind=IGA_INTEGER_KIND)  :: q
-  real   (kind=IGA_REAL_KIND   )  :: DetF
   do q=1,nqp
      call GeometryMap(&
           order,&
           nen,X,&
           M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
           N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          DetF,F(:,:,q),G(:,:,q))
-     DetJac(q) = DetJac(q) * DetF
+          DetF(q),F(:,:,q),G(:,:,q))
   end do
 contains
 include 'petigageo.f90.in'
   if (point->geometry) {
     PetscReal **M = point->basis;
     PetscReal **N = point->shape;
-    PetscReal dX  = point->detJac;
+    PetscReal *dX = &point->detX;
     PetscReal **gX = point->gradX;
     switch (point->dim) {
     case 3: IGA_ShapeFuns_3D(order,
                              point->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]); break;
     case 2: IGA_ShapeFuns_2D(order,
                              1,point->nen,
                              point->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]); break;
     case 1: IGA_ShapeFuns_1D(order,
                              1,point->nen,
                              point->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]); break;
     }
   }
   PetscFunctionReturn(0);
   ierr = PetscFree(element->basis[2]);CHKERRQ(ierr);
   ierr = PetscFree(element->basis[3]);CHKERRQ(ierr);
 
+  ierr = PetscFree(element->detX);CHKERRQ(ierr);
   ierr = PetscFree(element->gradX[0]);CHKERRQ(ierr);
   ierr = PetscFree(element->gradX[1]);CHKERRQ(ierr);
   ierr = PetscFree(element->shape[0]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen*dim*dim,PetscReal,&element->basis[2]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen*dim*dim*dim,PetscReal,&element->basis[3]);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*nen,PetscReal,&element->shape[0]);CHKERRQ(ierr);
     ierr = PetscMemzero(element->basis[2],sizeof(PetscReal)*nqp*nen*dim*dim);CHKERRQ(ierr);
     ierr = PetscMemzero(element->basis[3],sizeof(PetscReal)*nqp*nen*dim*dim*dim);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->shape[0],sizeof(PetscReal)*nqp*nen);CHKERRQ(ierr);
   }
   if (element->dim == element->nsd) /* XXX */
   if (element->geometry) {
+    PetscInt q, nqp = element->nqp;
     PetscReal **M = element->basis;
     PetscReal **N = element->shape;
-    PetscReal *dX = element->detJac;
+    PetscReal *dX = element->detX;
     PetscReal **gX = element->gradX;
     switch (element->dim) {
     case 3: IGA_ShapeFuns_3D(order,
                              N[0],N[1],N[2],N[3],
                              dX,gX[0],gX[1]); break;
     }
+    for (q=0; q<nqp; q++)
+      element->detJac[q] *= dX[q];
   }
   PetscFunctionReturn(0);
 }

src/petigaftn.F90

      type(C_PTR) :: point
      type(C_PTR) :: scale
      type(C_PTR) :: basis(0:3)
+     type(C_PTR) :: detX
      type(C_PTR) :: gradX(0:1)
      type(C_PTR) :: shape(0:3)
   end type IGAPoint

src/petigapoint.c

   point->basis[2] += nen*dim*dim;
   point->basis[3] += nen*dim*dim*dim;
 
+  point->detX     += 1;
   point->gradX[0] += dim*dim;
   point->gradX[1] += dim*dim;
   point->shape[0] += nen;
   point->basis[3] = element->basis[3];
 
   if (element->geometry && dim == nsd) { /* XXX */
+    point->detX     = element->detX;
     point->gradX[0] = element->gradX[0];
     point->gradX[1] = element->gradX[1];
     point->shape[0] = element->shape[0];
 #define __FUNCT__ "IGAPointAddMat"
 PetscErrorCode IGAPointAddMat(IGAPoint point,const PetscScalar k[],PetscScalar K[])
 {
-  PetscInt       n;
+  PetscInt       m,n;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(point,1);
   PetscValidScalarPointer(k,2);
   PetscValidScalarPointer(K,3);
-  n = point->nen*point->dof;
-  ierr = IGAPointAddArray(point,n*n,k,K);CHKERRQ(ierr);
+  m = n = point->nen*point->dof;
+  ierr = IGAPointAddArray(point,m*n,k,K);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
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.