Skip to main content

Assemble

This procedure pointer assembles the problem

Interface​

INTERFACE
  MODULE SUBROUTINE Assemble(obj)
    CLASS(AbstractSteadyStokes_), INTENT(INOUT) :: obj
  END SUBROUTINE Assemble
END INTERFACE

Example​

This example tests and demonstrates the usage of AssembleTanMat() method.

Use modules

PROGRAM main
  USE easifemBase
  USE easifemClasses
  USE easifemMaterials
  USE easifemKernels
  USE SteadyStokes221_Class
  IMPLICIT NONE

Declare variables

  TYPE( SteadyStokes221_ ) :: obj
  TYPE( ParameterList_ ) :: param
  TYPE( DomainPointer_ ) :: domains(2)
  CLASS( Domain_ ), POINTER :: dom => NULL()
  CHARACTER( LEN = * ), PARAMETER :: domainFileNamePressure="./mesh_tri3.h5"
  CHARACTER( LEN = * ), PARAMETER :: domainFileNameVelocity="./mesh_tri6.h5"
  TYPE(String) :: filename(2)
  TYPE( MeshSelection_ ) :: region
  CLASS( DirichletBC_ ), POINTER :: dbc => NULL()

Set parameters

Initiate an instance of ParameterList_, param, this will be used to initiate several objects.

Initiate param
  CALL FPL_INIT(); CALL param%Initiate()
set SteadyStokes221 parameters
  CALL SetSteadyStokes221Param( &
    & param=param, &
    & isConservativeForm=.TRUE., &
    & gravity = [0.0_DFP, -9.8_DFP, 0.0_DFP], &
    & domainFileForPressure = domainFileNamePressure, &
    & domainFileForVelocity = domainFileNameVelocity, &
    & engine="NATIVE_SERIAL", &
    & CoordinateSystem=KERNEL_CARTESIAN, &
    & maxIter = 100, &
    & rtoleranceForPressure = REAL( 1.0E-6, DFP ), &
    & rtoleranceForVelocity = REAL( 1.0E-6, DFP ), &
    & atoleranceForPressure = REAL( 1.0E-6, DFP ), &
    & atoleranceForVelocity = REAL( 1.0E-6, DFP ), &
    & toleranceForSteadyState = REAL( 1.0E-6, DFP ), &
    & tFluidMaterials=1, &
    & tDirichletBCForPressure=1, &
    & tDirichletBCForVelocity=3, &
    & baseInterpolationForSpace="LagrangeInterpolation", &
    & baseContinuityForSpace="H1", &
    & quadratureTypeForSpace="GaussLegendre", &
    & postProcessOpt=1)
Set param for linSolver
  CALL SetLinSolverParam( &
    & param=param, &
    & solverName=LIS_GMRES,&
    & preconditionOption=NO_PRECONDITION, &
    & convergenceIn=convergenceInRes, &
    & convergenceType=relativeConvergence, &
    & maxIter=100, &
    & relativeToRHS=.TRUE., &
    & KrylovSubspaceSize=20, &
    & rtol=REAL( 1.0E-10, DFP ), &
    & atol=REAL( 1.0D-10, DFP ) )

Initiate domain

Initiate domain
  filename = [String(domainFileNameVelocity), &
    & String(domainFileNamePressure)]
  CALL e%setQuietMode(.TRUE.)
  CALL Initiate(domains=domains, filename=filename)
In line 1, make sure that Velocity domain comes first.

Initiate kernel

Initiate kernel
  CALL e%setQuietMode(.TRUE.)
  CALL obj%Initiate(param=param, domains=domains )

Here, line 1 is to allow verbosity.

Add fluid material to kernel. To do so, we first create an instance of MeshSelection_. Then we add this instance to the kernel.

adding fluid material
  CALL e%setQuietMode(.TRUE.)
  CALL region%Initiate( isSelectionByMeshID=.TRUE. )
  CALL region%Add( dim=obj%nsd, meshID=[1] )
  CALL SetFluidMaterialParam( &
    & param=param, &
    & name="fluidMaterial", &
    & massDensity=1000.0_DFP, &
    & dynamicViscosity = 0.001_DFP, &
    & stressStrainModel="NewtonianFluidModel" )
  CALL SetNewtonianFluidModelParam( &
    & param = param, &
    & dynamicViscosity = 0.001_DFP )
  CALL e%setQuietMode(.TRUE.)
  CALL obj%AddFluidMaterial( &
    & materialNo=1, &
    & materialName="fluidMaterial", &
    & param=param, &
    & region=region)
  CALL region%Deallocate()

Add Dirichlet boundary condition

V1=0, Now we show how to add dirichlet boundary condition. To this end first we create an instance of MeshSelection_ to select the region of the mesh. Then we define the dirichlet bonundary condition, and pass these two information to kernel.

Set parameters for dirichlet boundary condition:

  CALL SetDirichletBCParam( &
    & param=param, &
    & name="ZeroV1", &
    & idof=1, &
    & nodalValueType=Constant, &
    & useFunction=.FALSE. )

Select the mesh region:

#define BOTTOM 1
#define RIGHT 2
#define TOP 3
#define LEFT 4
  CALL region%Initiate( isSelectionByMeshID=.TRUE. )
  CALL region%Add( dim=obj%nsd-1, meshID=[BOTTOM, RIGHT, LEFT] )
  CALL region%Set()

Add dirichlet boundary condition and the region to kernel:

  CALL obj%AddVelocityDirichletBC( &
    & dbcNo=1, &
    & param=param, &
    & boundary=region )
  dbc => obj%GetVelocityDirichletBCPointer( dbcNo=1 )
  CALL dbc%Set( ConstantNodalValue=0.0_DFP )
  dbc=>NULL()

V1=U, Set parameters for dirichlet boundary condition:

  CALL SetDirichletBCParam( &
    & param=param, &
    & name="UpstreamV1", &
    & idof=1, &
    & nodalValueType=Constant, &
    & useFunction=.FALSE. )

Select the mesh region:

  CALL region%Deallocate()
  CALL region%Initiate( isSelectionByMeshID=.TRUE. )
  CALL region%Add( dim=obj%nsd-1, meshID=[TOP] )
  CALL region%Set()

Add dirichlet boundary condition and the region to kernel:

  CALL obj%AddVelocityDirichletBC( &
    & dbcNo=2, &
    & param=param, &
    & boundary=region )
  dbc => obj%GetVelocityDirichletBCPointer( dbcNo=2 )
  CALL dbc%Set( ConstantNodalValue=0.1_DFP )
  dbc=>NULL()
  CALL region%Deallocate()

V2=0, Set parameters for dirichlet boundary condition:

  CALL SetDirichletBCParam( &
    & param=param, &
    & name="ZeroV2", &
    & idof=2, &
    & nodalValueType=Constant, &
    & useFunction=.FALSE. )

Select the mesh region:

  CALL region%Deallocate()
  CALL region%Initiate( isSelectionByMeshID=.TRUE. )
  CALL region%Add( dim=obj%nsd-1, meshID=[BOTTOM, RIGHT, TOP, LEFT] )
  CALL region%Set()

Add dirichlet boundary condition and the region to kernel:

  CALL obj%AddVelocityDirichletBC( &
    & dbcNo=3, &
    & param=param, &
    & boundary=region )
  dbc => obj%GetVelocityDirichletBCPointer( dbcNo=3 )
  CALL dbc%Set( ConstantNodalValue=0.0_DFP ); dbc=>NULL()
  CALL region%Deallocate()

P=0, Set parameters for dirichlet boundary condition:

setting boundary condition P=0
CALL SetDirichletBCParam( &
  & param=param, &
  & name="ZeroP", &
  & idof=1, &
  & nodalValueType=Constant, &
  & useFunction=.FALSE. )

Select the mesh region:

CALL region%Deallocate()
CALL region%Initiate( isSelectionByMeshID=.TRUE. )
CALL region%Add( dim=obj%nsd-1, meshID=[BOTTOM] )
CALL region%Set()

Add dirichlet boundary condition and the region to kernel:

  CALL obj%AddPressureDirichletBC( dbcNo=1, param=param, &
    & boundary=region )
  dbc => obj%GetPressureDirichletBCPointer( dbcNo=1 )
  CALL dbc%Set( ConstantNodalValue=0.0_DFP ); dbc=>NULL()
  CALL region%Deallocate()

Now that we are done with the setup, we should call Set method. In this method, the kernel checks the data, configuration, and intiates appropriate variables.

  CALL obj%Set()

AssembleTanMat

assembling tangent matrix
  CALL e%setQuietMode(.TRUE.)
  CALL obj%AssembleTanmat()
checking sparsity pattern
  CALL obj%Kmat%SPY(filename="./Kmat", ext=".png")
  CALL obj%Gmat%SPY(filename="./Gmat", ext=".png")

Results are given below

Kmat

Assemble

assemble
  CALL e%setQuietMode(.FALSE.)
  CALL obj%Assemble()

Cleanup

  CALL obj%Deallocate()
  CALL param%Deallocate()
  CALL FPL_FINALIZE()
END PROGRAM main