EbOM_example_mech_3bars.f90 Source File

program main
use SparseClass
use FEMDomainClass
use FEMSolverClass
implicit none


type(FEMDomain_),target :: cube(3)
type(FEMDomainp_) :: cube_p(3)
type(FEMSolver_) :: solver
integer(int32),allocatable   :: FixBoundary(:)
integer(int32),allocatable   :: Overset_Elements(:,:)
integer(int32)   :: ElementID, DomainID,i
type(IO_) :: f
type(MPI_) :: mpid

integer(int32),allocatable :: row(:)
integer(int32),allocatable :: col(:)
real(real64),allocatable :: val(:),all_disp(:)
integer(int32) :: nn = 1000000.0

call mpid%start()

! mesh
call cube(1)%create("Cube3D",x_num=30,y_num=3,z_num=3)
call cube(2)%create("Cube3D",x_num=30,y_num=3,z_num=3)
call cube(3)%create("Cube3D",x_num=30,y_num=3,z_num=3)

call cube(1)%resize(x=10.0d0)
call cube(2)%resize(x=10.0d0)
call cube(3)%resize(x=10.0d0)

call cube(2)%move(x=cube(1)%x_max()-0.10d0 )
call cube(3)%move(x=cube(2)%x_max()-0.10d0 )

call cube(1)%vtk("cube_no1")
call cube(2)%vtk("cube_no2")
call cube(3)%vtk("cube_no3")

print *, "[ok] Mesh creation"

! overset
call cube(2)%overset(cube(3),DomainID=3, &
    algorithm=FEMDomain_Overset_GPP,MyDomainID=2) ! or "P2P"

call cube(2)%overset(cube(1),DomainID=1, &
    algorithm=FEMDomain_Overset_GPP,MyDomainID=2) ! or "P2P"

print *, cube(2)%numOversetElements()

cube_p(1)%femdomainp => cube(1)
cube_p(2)%femdomainp => cube(2)
cube_p(3)%femdomainp => cube(3)

call solver%init(NumDomain=3)
call solver%setDomain(FEMDomainPointers=cube_p(:))
call solver%setCRS(DOF=3)

call print(solver%CRS_ID_Starts_From ) 

!$OMP parallel 
!$OMP do
do DomainID=1,3
    do ElementID = 1, cube(DomainID)%ne()
        call solver%setMatrix(DomainID=DomainID,ElementID=ElementID,DOF=3,&
           Matrix=cube(DomainID)%StiffnessMatrix(ElementID=ElementID,E=5000000.0d0, v=0.330d0) )
        call solver%setVector(DomainID=DomainID,ElementID=ElementID,DOF=3,&
            Vector=cube(DomainID)%MassVector(&
                ElementID=ElementID,&
                DOF=cube(DomainID)%nd() ,&
                Density=1.700d0,&
                Accel=[0.0d0, 0.0d0, -9.80d0]&
                ) & 
            )
    enddo
enddo
!$OMP end do
!$OMP end parallel

print *, "[ok]Element-matrices done"

call solver%setEbOM(penalty=10000000.0d0, DOF=3)

!print *, "matrices imported."
! disp. boundary

FixBoundary = cube(1)%select(x_max = cube(1)%x_min() )*3-2
call solver%fix(DomainID=1,IDs=FixBoundary,FixValue=0.0d0)
FixBoundary(:) = FixBoundary(:) +1
call solver%fix(DomainID=1,IDs=FixBoundary,FixValue=0.0d0)
FixBoundary(:) = FixBoundary(:) +1
call solver%fix(DomainID=1,IDs=FixBoundary,FixValue=0.0d0)
!FixBoundary = cube(3)%select(x_min = cube(3)%x_max() )*3-0
!FixBoundary(:) = FixBoundary(:) + solver%CRS_ID_Starts_From(3)-1
!call solver%fix(DomainID=3,IDs=FixBoundary,FixValue=-3.00d0)
!
print *, "b.c. imported."
!
call cube(1)%vtk("cube_1_step_0")
call cube(2)%vtk("cube_2_step_0")
call cube(3)%vtk("cube_3_step_0")
!
!! solve
solver%debug = .true.
!
!solver%er0 = 
solver%itrmax = 10000
all_disp = solver%solve()
call cube(1)%deform(disp=all_disp(1:solver%CRS_ID_Starts_From(2)-1  ) )
call cube(2)%deform(disp=all_disp(solver%CRS_ID_Starts_From(2):solver%CRS_ID_Starts_From(3)-1  ) )
call cube(3)%deform(disp=all_disp(solver%CRS_ID_Starts_From(3):  ) )
!
!
call cube(1)%vtk("cube_1_step_1" )
call cube(2)%vtk("cube_2_step_1" )
call cube(3)%vtk("cube_3_step_1" )
!
call mpid%end()


end program