3.1 IFEM_Solve_shared

void IFEM_Solve_shared(ILSI_Solver s,ILSI_Param *p, int fem_mesh, int fem_entity,int length,int width, IFEM_Matrix_product_c A, void *ptr, const double *b, double *x);
subroutine IFEM_Solve_shared(s,p, fem_mesh,fem_entity,length,width, A,ptr,b,x)
external solver subroutine :: s
double precision, intent(inout) :: p(ILSI_PARAM)
integer, intent(in) :: fem_mesh, fem_entity, length,width
external matrix-vector product subroutine :: A
TYPE(varies), pointer :: ptr
double precision, intent(in) :: b(width,length)
double precision, intent(inout) :: x(width,length)

This routine solves the linear system for the unknown vector . s and p give the particular linear solver to use, and are described in more detail in Section 2. fem_mesh and fem_entity give the FEM framework mesh (often FEM_Mesh_default_read()) and entity (often FEM_NODE) with which the known and unknown vectors are listed.

width gives the number of degrees of freedom (entries in the vector) per node. For example, if there is one degree of freedom per node, width is one. length should always equal the number of FEM nodes.

A is a local matrix-vector product routine you must write. Its interface is described in Section 3.1.1. ptr is a pointer passed down to A-it is not otherwise used by the framework.

b is the known vector. x, on input, is the initial guess for the unknown vector. On output, x is the final value for the unknown vector. b and x should both have length * width entries. In C, DOF of node should be indexed as width. In Fortran, these arrays should be allocated like x(width,length).

When this routine returns, x is the final value for the unknown vector, and the output values of the solver parameters p will have been written.

// C++ Example
  int mesh=FEM_Mesh_default_read();
  int nNodes=FEM_Mesh_get_length(mesh,FEM_NODE);
  int width=3; //A 3D problem
  ILSI_Param solverParam;
  struct myProblemData myData;
  
  double *b=new double[nNodes*width];
  double *x=new double[nNodes*width];
  ... prepare solution target b and guess x ...
  
  ILSI_Param_new(&solverParam);
  solverParam.maxResidual=1.0e-4;
  solverParam.maxIterations=500;
  
  IFEM_Solve_shared(IFEM_CG_Solver,&solverParam,
         mesh,FEM_NODE, nNodes,width,
         myMatrixVectorProduct, &myData, b,x);
  
! F90 Example
  include 'ifemf.h'
  INTEGER :: mesh, nNodes,width
  DOUBLE PRECISION, ALLOCATABLE :: b(:,:), x(:,:)
  DOUBLE PRECISION :: solverParam(ILSI_PARAM)
  TYPE(myProblemData) :: myData
  
  mesh=FEM_Mesh_default_read()
  nNodes=FEM_Mesh_get_length(mesh,FEM_NODE)
  width=3   ! A 3D problem
  
  ALLOCATE(b(width,nNodes), x(width,nNodes))
  ... prepare solution target b and guess x ..
  
  ILSI_Param_new(&solverParam);
  solverParam(1)=1.0e-4;
  solverParam(2)=500;
  
  IFEM_Solve_shared(IFEM_CG_Solver,solverParam,
         mesh,FEM_NODE, nNodes,width,
         myMatrixVectorProduct, myData, b,x);

3.1.1 Matrix-vector product routine

IFEM requires you to write a matrix-vector product routine that will evaluate for various vectors . You may use any subroutine name, but it must take these arguments:

void IFEM_Matrix_product(void *ptr,int length,int width, const double *src, double *dest);
subroutine IFEM_Matrix_product(ptr,length,width,src,dest)
TYPE(varies), pointer :: ptr
integer, intent(in) :: length,width
double precision, intent(in) :: src(width,length)
double precision, intent(out) :: dest(width,length)

The framework calls this user-written routine when it requires a matrix-vector product. This routine should compute , interpreting and as vectors. length gives the number of nodes and width gives the number of degrees of freedom per node, as above.

In writing this routine, you are responsible for choosing a representation for the matrix . For many problems, there is no need to represent explicitly-instead, you simply evaluate by looping over local elements, taking into account the values of . This example shows how to write the matrix-vector product routine for simple 1D linear elastic springs, while solving for displacement given net forces.

After calling this routine, the framework will handle combining the overlapping portions of these vectors across processors to arrive at a consistent global matrix-vector product.

// C++ Example
#include "ifemc.h"

typedef struct {
  int nElements; //Number of local elements
  int *conn; // Nodes adjacent to each element: 2*nElements entries
  double k; //Uniform spring constant
} myProblemData;

void myMatrixVectorProduct(void *ptr,int nNodes,int dofPerNode,
          const double *src,double *dest)
{
  myProblemData *d=(myProblemData *)ptr;
  int n,e;
  // Zero out output force vector:
  for (n=0;n<nNodes;n++) dest[n]=0;
  // Add in forces from local elements
  for (e=0;e<d->nElements;e++)
    int n1=d->conn[2*e+0]; // Left node
    int n2=d->conn[2*e+1]; // Right node
    double f=d->k * (src[n2]-src[n1]); //Force
    dest[n1]+=f;
    dest[n2]-=f;
  
}


! F90 Example
  TYPE(myProblemData)
    INTEGER :: nElements
    INTEGER, ALLOCATABLE :: conn(2,:)
    DOUBLE PRECISION :: k
  END TYPE
  
SUBROUTINE myMatrixVectorProduct(d,nNodes,dofPerNode,src,dest)
  include 'ifemf.h'
  TYPE(myProblemData), pointer :: d
  INTEGER :: nNodes,dofPerNode
  DOUBLE PRECISION :: src(dofPerNode,nNodes), dest(dofPerNode,nNodes)
  INTEGER :: e,n1,n2
  DOUBLE PRECISION :: f
  
  dest(:,:)=0.0
  do e=1,d%nElements
    n1=d%conn(1,e)
    n2=d%conn(2,e)
    f=d%k * (src(1,n2)-src(1,n1))
    dest(1,n1)=dest(1,n1)+f
    dest(1,n2)=dest(1,n2)+f
  end do
END SUBROUTINE