Updated for new FEM implementation.
authorOrion Lawlor <olawlor@acm.org>
Mon, 6 Nov 2000 03:38:35 +0000 (03:38 +0000)
committerOrion Lawlor <olawlor@acm.org>
Mon, 6 Nov 2000 03:38:35 +0000 (03:38 +0000)
tests/fem/femtest/cmesh.dat
tests/fem/femtest/fpgm.f90
tests/fem/femtest/pgm.C

index e39f1dcb2830b5eaf3d7c97cc9bd9d013ec77c0f..68763324f94c85acf91eacb53a14f212ec749508 100644 (file)
@@ -1,3 +1,15 @@
+9 16 4
+0 1 4 5
+1 2 5 6
+2 3 6 7
+4 5 8 9
+5 6 9 10
+6 7 10 11
+8 9 12 13
+9 10 13 14 
+10 11 14 15
+
+# 3x3 mesh:
 4 9 4
 0 1 3 4
 1 2 4 5
index 17219b69986d551b32757c34977ef6da5fb5e1cf..b8bd6f42ef5806a0c789412b1d5a20e2218a3968 100644 (file)
@@ -2,50 +2,66 @@ subroutine init()
 implicit none
 include 'femf.h'
 
-  integer :: i, j, nelems, nnodes, ctype, esize
+  integer :: i, j, nelems, nnodes, esize
   integer, dimension(:,:), allocatable:: conn
+  double precision, dimension(:,:), allocatable:: nodeData
 
   call FEM_Print('init called')
   open(20, file='fmesh.dat')
-  read(20,*) nelems, nnodes, ctype
-  if (ctype .eq. FEM_TRIANGULAR) then
-    esize = 3
-  else 
-    if(ctype .eq. FEM_HEXAHEDRAL) then
-      esize = 8
-    else
-      esize = 4
-    endif
-  endif
+  read(20,*) nelems, nnodes, esize
   allocate(conn(nelems, esize))
   do i=1,nelems
     read(20,*) (conn(i,j),j=1,esize)
   enddo
   close(20)
-  call FEM_Set_Mesh(nelems, nnodes, ctype, conn)
+  call FEM_Set_Elem(1,nelems,0,esize)
+  call FEM_Set_Elem_Conn_c(1,conn)
+  
+  allocate(nodeData(nnodes,2))
+  do i=1,nnodes
+     nodeData(i,1)=0
+     nodeData(i,2)=0
+  enddo
+  nodeData(1,1)=1
+  nodeData(1,2)=0.25
+  nodeData(2,2)=0.25
+  nodeData(4,2)=0.25
+  nodeData(5,2)=0.25
+  call FEM_Set_Node(nnodes,2)
+  call FEM_Set_Node_Data_c(nodeData)
 end subroutine init
 
-subroutine driver(nnodes, nnums, nelems, enums, npere, conn)
+subroutine driver()
 implicit none
 include 'femf.h'
 
-  integer  :: nnodes, nelems, npere
-  integer, dimension(nnodes) :: nnums
-  integer, dimension(nelems) :: enums
-  integer, dimension(nelems, npere) :: conn
-
+  integer  :: nnodes,nodeDataP, nelems, elemData,npere
+  integer, dimension(:,:), allocatable:: conn
+  double precision, dimension(:,:), allocatable:: nodeData
   integer :: i, j, fid
   logical :: failed
   double precision :: sum
-  double precision, dimension(nnodes) :: nodes
-  double precision, dimension(nelems) :: elements
+  double precision, dimension(:), allocatable:: nodes
+  double precision, dimension(:), allocatable:: elements
+
 
-  !call FEM_Print_Partition()
+  call FEM_Get_Elem(1,nelems,elemData,npere)
+  call FEM_Get_Node(nnodes,nodeDataP)
+
+  allocate(conn(nelems, npere))
+  call FEM_Get_Elem_Conn_c(1,conn)
+  allocate(nodeData(nnodes,nodeDataP))
+  call FEM_Get_Node_Data_c(nodeData);
+
+  allocate(nodes(nnodes))
+  allocate(elements(nelems))
+
+  call FEM_Print_Partition()
 
   nodes = 0.0
   elements = 0.0
   do i=1,nnodes
-    if (nnums(i) .eq. 1) nodes(i) = 1.0
+     nodes(i)=nodeData(i,1)
   enddo
   fid = FEM_Create_Field(FEM_DOUBLE, 1, 0, 8)
   do i=1,nelems
@@ -63,12 +79,7 @@ include 'femf.h'
   call FEM_Update_Field(fid, nodes(1))
   failed = .FALSE.
   do i=1,nnodes
-    if (nnums(i).eq.1 .or. nnums(i).eq.2 .or. &
-&       nnums(i).eq.4 .or. nnums(i).eq.5) then 
-      if(nodes(i) .ne. 0.25) failed = .TRUE.
-    else
-      if (nodes(i) .ne. 0.0) failed = .TRUE.
-    endif
+    if (nodes(i) .ne. nodeData(i,2)) failed= .TRUE.
   enddo
   if (failed) then
     call FEM_Print('update_field test failed.')
@@ -92,6 +103,13 @@ include 'femf.h'
   call FEM_Done()
 end subroutine driver
 
+subroutine mesh_updated(param)
+  implicit none
+  integer :: param
+  include 'femf.h'
+  call FEM_Print('mesh_updated called')
+end subroutine
+
 subroutine finalize()
 implicit none
 include 'femf.h'
index 3f9bc66643bda29343aa7cbb950cfda0d41be7fe..9e0e0e0196ded7ff3e9d6181b1b968af7a0a6272 100644 (file)
@@ -1,28 +1,80 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "charm++.h"
 #include "fem.h"
 
+//These values are the sum of the node values
+// after each timestep.
+//If you change the mesh size, initial conditions, etc.
+//  you'll have to update this array with the correct values,
+//  which are (conveniently) printed out by init.
+const double reduceValues[4]={512.00, 512.00, 465.00, -1};
+
 extern "C" void
 init(void)
 {
   CkPrintf("init called\n");
-  FILE *fp = fopen("cmesh.dat", "r");
-  if(fp==0) { CkAbort("Cannot open cmesh.dat file.\n"); }
-  int nelems, nnodes, ctype;
-  fscanf(fp, "%d%d%d\n", &nelems, &nnodes, &ctype);
-  int i,j;
-  int esize = (ctype==FEM_TRIANGULAR)?3:
-              ((ctype==FEM_HEXAHEDRAL)?8:
-             4);
-  int *conn = new int[nelems*esize];
-  for(i=0;i<nelems;i++) {
-    for(j=0;j<esize;j++) {
-      fscanf(fp,"%d",&conn[i*esize+j]);
+  const int tsteps=3;//Number of time steps to simulate
+  const int dim=5;//Length of one side of the FEM mesh
+  const int np=4; //Nodes per element
+  int conn[dim*dim*np];
+  double elements[dim*dim]={0.0};
+  double nodes[(dim+1)*(dim+1)]={0.0};
+  double noData[(dim+1)*(dim+1)*tsteps];
+
+  int nelems=dim*dim, nnodes=(dim+1)*(dim+1);
+
+  //Describe the nodes and elements
+  FEM_Set_Node(nnodes,tsteps);
+  FEM_Set_Elem(0,nelems,1,np);
+
+  //Create the connectivity array
+  for(int y=0;y<dim;y++) for (int x=0;x<dim;x++) {
+          conn[(y*dim+x)*np+0]=(y  )*(dim+1)+(x  );
+          conn[(y*dim+x)*np+1]=(y+1)*(dim+1)+(x  );
+          conn[(y*dim+x)*np+2]=(y+1)*(dim+1)+(x+1);
+          conn[(y*dim+x)*np+3]=(y  )*(dim+1)+(x+1);
+  }
+
+  //Set the initial conditions
+  elements[3*dim+1]=256;
+  elements[2*dim+1]=256;
+  FEM_Set_Elem_Data(0,elements);
+  FEM_Set_Elem_Conn(0,conn);
+
+  //Run the time loop over our serial mesh--
+  // we'll use this data to check the parallel calculation.
+  CkPrintf("const double reduceValues[%d]={",tsteps+1);
+  for (int t=0;t<tsteps;t++)
+  {
+    int i,j;
+
+       //Nodes are sum of surrounding elements
+    for(i=0;i<nnodes;i++) nodes[i] = 0.0;
+    for(i=0;i<nelems;i++)
+      for(j=0;j<np;j++)
+        nodes[conn[i*np+j]] += elements[i]/np;
+
+       //Elements are average of surrounding nodes
+    for(i=0;i<nelems;i++) {
+         double sum=0;
+      for(j=0;j<np;j++)
+        sum += nodes[conn[i*np+j]];
+      elements[i] = sum/np;
     }
+       
+    //Save the node values for this timestep
+       for (i=0;i<nnodes;i++)
+                noData[i*tsteps+t]=nodes[i];      
+
+       //Compute the sum across all nodes
+       double reduceSum=0;
+       for (i=0;i<nnodes;i++) reduceSum+=nodes[i];
+    CkPrintf("%.2f, ",reduceSum);
   }
-  fclose(fp);
-  FEM_Set_Mesh(nelems, nnodes, ctype, conn);
+  FEM_Set_Node_Data(noData);
+  CkPrintf("-1};\n");
 }
 
 typedef struct _node {
@@ -34,63 +86,111 @@ typedef struct _element {
 } Element;
 
 extern "C" void
-driver(int nnodes, int *nnums, int nelems, int *enums, int npere, int *conn)
+driver(void)
 {
+  int nnodes,nelems,nnodeData,nelemData,np;
+
+  FEM_Get_Node(&nnodes,&nnodeData);
+  double *nodeData=new double[nnodeData*nnodes];
+  FEM_Get_Node_Data(nodeData);  
+
+  FEM_Get_Elem(0,&nelems,&nelemData,&np);
+  int *conn=new int[np*nelems];
+  FEM_Get_Elem_Conn(0,conn);
+  double *elData=new double[nelemData*nelems];
+  FEM_Get_Elem_Data(0,elData);
+
+  int tsteps=nnodeData;
+
   int myId = FEM_My_Partition();
-  // FEM_Print_Partition();
+  FEM_Print_Partition();
   Node *nodes = new Node[nnodes];
   Element *elements = new Element[nelems];
   int i;
+
+//Set initial conditions
   for(i=0;i<nnodes;i++) {
-    nodes[i].val = 0.0;
-    if(nnums[i]==0) { nodes[i].val = 1.0; }
+    nodes[i].val = 0;
   }
-  for(i=0;i<nelems;i++) { elements[i].val = 0.0; }
+  for(i=0;i<nelems;i++) { elements[i].val = elData[i]; }
   int fid = FEM_Create_Field(FEM_DOUBLE, 1, 0, sizeof(Node));
-  int j;
-  for(i=0;i<nelems;i++) {
-    for(j=0;j<npere;j++) {
-      elements[i].val += nodes[conn[i*npere+j]].val;
+
+//Time loop
+  for (int t=0;t<tsteps;t++)
+  {
+    int i,j;
+
+       //Nodes are sum of surrounding elements
+    for(i=0;i<nnodes;i++) nodes[i].val = 0.0;
+    for(i=0;i<nelems;i++)
+      for(j=0;j<np;j++)
+        nodes[conn[i*np+j]].val += elements[i].val/np;
+
+       //Update shared nodes
+    FEM_Update_Field(fid, nodes);
+
+       //Elements are average of surrounding nodes
+    for(i=0;i<nelems;i++) {
+         double sum=0;
+      for(j=0;j<np;j++)
+        sum += nodes[conn[i*np+j]].val;
+      elements[i].val = sum/np;
     }
-    elements[i].val /= npere;
-  }
-  for(i=0;i<nnodes;i++) { nodes[i].val = 0.0; }
-  for(i=0;i<nelems;i++) {
-    for(j=0;j<npere;j++) {
-      nodes[conn[i*npere+j]].val += elements[i].val;
+
+       //Check the update
+    int failed = 0;
+    for(i=0;i<nnodes;i++)
+        if(nodes[i].val != nodeData[nnodeData*i+t]) 
+                        failed = 1;
+    
+    if(failed==0) {
+      CkPrintf("[chunk %d] update_field %d test passed.\n", myId, t);
+    } else {
+      CkPrintf("[chunk %d] update_field %d test failed.\n", myId, t);
     }
-  }
-  FEM_Update_Field(fid, nodes);
-  int failed = 0;
-  for(i=0;i<nnodes;i++) {
-    if(nnums[i]==0 || nnums[i]==1 || nnums[i]==3 || nnums[i]==4) {
-      if(nodes[i].val != 0.25) { failed = 1; }
+
+    double sum = 0.0;
+    FEM_Reduce_Field(fid, nodes, &sum, FEM_SUM);
+    if(sum==reduceValues[t]) {
+      CkPrintf("[chunk %d] reduce_field test passed.\n", myId);
     } else {
-      if(nodes[i].val != 0.0) { failed = 1; }
+      CkPrintf("[chunk %d] reduce_field test failed (sum=%f at %d).\n", myId, sum, t);
     }
+
+    FEM_Set_Node(nnodes,1);
+    FEM_Set_Node_Data((double *)nodes);
+    FEM_Update_Mesh(1+t,0);
   }
-  if(failed==0) {
-    CkPrintf("[chunk %d] update_field test passed.\n", myId);
-  } else {
-    CkPrintf("[chunk %d] update_field test failed.\n", myId);
-  }
-  double sum = 0.0;
-  FEM_Reduce_Field(fid, nodes, &sum, FEM_SUM);
-  if(sum==1.0) {
-    CkPrintf("[chunk %d] reduce_field test passed.\n", myId);
-  } else {
-    CkPrintf("[chunk %d] reduce_field test failed.\n", myId);
-  }
-  sum = 1.0;
-  FEM_Reduce(fid, &sum, &sum, FEM_SUM);
-  if(sum==(double)FEM_Num_Partitions()) {
+
+  double localSum = 1.0,globalSum;
+  FEM_Reduce(fid, &localSum, &globalSum, FEM_SUM);
+  if(globalSum==(double)FEM_Num_Partitions()) {
     CkPrintf("[chunk %d] reduce test passed.\n", myId);
   } else {
     CkPrintf("[chunk %d] reduce test failed.\n", myId);
   }
+
   FEM_Done();
 }
 
+extern "C" void
+mesh_updated(int param)
+{
+  int nnodes,dataPer;
+  CkPrintf("mesh_updated(%d) called.\n",param);
+  FEM_Get_Node(&nnodes,&dataPer);
+  double *ndata=new double[nnodes*dataPer];
+  FEM_Get_Node_Data(ndata);
+  double sum=0;
+  for (int i=0;i<nnodes;i++)
+    sum+=ndata[i];
+  if (sum!=reduceValues[param-1])
+    CkPrintf("mesh_updated test FAILED! expected %f; got %f!\n",reduceValues[param-1],sum);
+  else
+    CkPrintf("mesh_updated test passed for %d\n",param);
+  delete[] ndata;
+}
+
 extern "C" void
 finalize(void)
 {