Update module dependencies

app/main.f90 → app/FE_test_sparsematrix.f90

@@ -1,7 +1,5 @@
 program main
-
-   use linear_pack, only: general_square_matrix
-   use quadrature_module, only: wp => quadrature_wp
+   use mod_kinds
    use M_attr
    use mod_FE_2D
 
@@ -37,105 +35,112 @@ end function cofunc_2d_t
    !> input information of field
    type(field) :: field_info
 
-   !> Field parameter input.
+!> Field parameter input.
    !>>TODO: read from toml file.
    real(wp) :: left = 0.0_wp, right = 10.0_wp, bottom = 0.0_wp, top = 10.0_wp
-      !! The problem domain is [left,right]*[bottom,top].
+   !! The problem domain is [left,right]*[bottom,top].
    integer :: Nh_parition = 100, Nv_parition = 100
-      !! The number of partition in horizontal and vertical direction.
+   !! The number of partition in horizontal and vertical direction.
    ! real(wp) :: h_partition, v_partition
-      !! The step size of the partition.
+   !! The step size of the partition.
    integer :: Nh_basis, Nv_basis
-      !! The number of FE basis functions in horizontal and vertical direction.
+   !! The number of FE basis functions in horizontal and vertical direction.
    integer :: Gauss_point_number = 6
-      !! The number of Gauss Quadrature points.
+   !! The number of Gauss Quadrature points.
    real(wp) :: tolerance = 1.0e-3_wp
-      !! The tolerance of the error.
+   !! The tolerance of the error.
    !> basis_type: the type of the FE.
    integer :: basis_type = 201
-      !! 201: 2D linear
-      !! 202: 2D quadratic
+   !! 201: 2D linear
+   !! 202: 2D quadratic
    integer :: mesh_type = 503
-      !! 503: 2D triangular mesh;
-      !! 504: 2D quadrilateral mesh; !!!NOT IMPLEMENTED YET!!!
+   !! 503: 2D triangular mesh;
+   !! 504: 2D quadrilateral mesh; !!!NOT IMPLEMENTED YET!!!
 
    real(wp), allocatable :: M(:, :), M_basis(:, :)
    integer, allocatable :: T(:, :), T_basis(:, :)
    integer, allocatable :: boundarynodes(:, :)
    integer, allocatable :: boundaryedges(:, :)
 
-   real(wp), allocatable :: A(:, :), b(:, :), A1(:, :), A2(:, :)
-   real(wp), allocatable :: solution(:, :)
+   real(wp), allocatable :: A(:, :), A1(:, :), A2(:, :)
+   real(wp), allocatable :: b(:), solution(:)
    ! real(wp) :: vet(2, 3)
 
+   !> sparse matrix and solver
+   type(sparse_COO):: sparse_A, sparseA1, sparseA2
+
+!> co-Function input.
    procedure(local_basis_func), pointer :: lbfunc => null()
    procedure(cofunc_2d), pointer :: cofunc => null()
    procedure(cofunc_2d), pointer :: cofuncf => null()
    lbfunc => local_basis_func_2D
    cofunc => function_c
    cofuncf => function_f
 
+!> Field information initialization.
    select case (basis_type)
-   case (201)
+    case (201)
       Nh_basis = Nh_parition
       Nv_basis = Nv_parition
-   case (202)
+    case (202)
       Nh_basis = 2*Nh_parition
       Nv_basis = 2*Nv_parition
    end select
 
    call field_info%init(left, right, bottom, top, &
-                        Nh_parition, Nv_parition, &
-                        Nh_basis, Nv_basis, &
-                        Gauss_point_number=Gauss_point_number, &
-                        tolerance=tolerance, &
-                        trial_basis_type=basis_type, &
-                        test_basis_type=basis_type, &
-                        mesh_type=mesh_type)
+      Nh_parition, Nv_parition, &
+      Nh_basis, Nv_basis, &
+      Gauss_point_number=Gauss_point_number, &
+      tolerance=tolerance, &
+      trial_basis_type=basis_type, &
+      test_basis_type=basis_type, &
+      mesh_type=mesh_type)
 
    ! call field_info%print('(A,2F6.3)')
 
    !> Allocate memory to A, b and solution. Set to zero.
    allocate (A(field_info%number_of_nodes_fe, field_info%number_of_nodes_fe))
    allocate (A1(field_info%number_of_nodes_fe, field_info%number_of_nodes_fe))
    allocate (A2(field_info%number_of_nodes_fe, field_info%number_of_nodes_fe))
-   allocate (b(field_info%number_of_nodes_fe, 1))
-   allocate (solution(field_info%number_of_nodes_fe, 1))
+   allocate (b(field_info%number_of_nodes_fe))
+   allocate (solution(field_info%number_of_nodes_fe))
    A1 = 0
    A2 = 0
    A = 0
    b = 0
    solution = 0
 
+   !> Generate the mesh and basis information matrix.
    call generate_info_matrix(field_info, M, T, &
-                             basis_type=basis_type, verbose=.true.)
+      basis_type=basis_type, verbose=.true.)
    call generate_info_matrix(field_info, M_basis, T_basis, &
-                             basis_type=basis_type, verbose=.true.)
+      basis_type=basis_type, verbose=.true.)
+   !> Generate the boundary nodes and edges matrix.
    call generate_boundarynodes(field_info, boundarynodes, verbose=.true.)
    call generate_boundaryedges(field_info, boundaryedges, verbose=.true.)
 
    !! Assemble the matrix and vector.
    !! Triangular mesh.
    call assemble_matirx_2D(A1, cofunc, lbfunc, &
-                           M, T, &
-                           T_basis, T_basis, &
-                           1, 0, &
-                           1, 0, &
-                           field_info)
+      M, T, &
+      T_basis, T_basis, &
+      1, 0, &
+      1, 0, &
+      field_info)
    call assemble_matirx_2D(A2, cofunc, lbfunc, &
-                           M, T, &
-                           T_basis, T_basis, &
-                           0, 1, &
-                           0, 1, &
-                           field_info)
+      M, T, &
+      T_basis, T_basis, &
+      0, 1, &
+      0, 1, &
+      field_info)
 
    A = A1 + A2
    deallocate (A1, A2)
    call assemble_vector_2D(b, cofuncf, lbfunc, &
-                           M, T, &
-                           T_basis, &
-                           0, 0, &
-                           field_info)
+      M, T, &
+      T_basis, &
+      0, 0, &
+      field_info)
    call treat_Dirchlet_boundary(function_g, A, b, boundarynodes, M)
 
    !! Solve the linear system.

example/example_1D/example_FE_1D.f90

@@ -1,6 +1,5 @@
 program example_FE_1D
-   use linear_pack, only: general_square_matrix
-   use quadrature_module, only: wp => quadrature_wp
+   use mod_kinds, only: wp
    use mod_FE_1D
 
    implicit none

example/example_2D/example_FE_2D.f90

@@ -1,6 +1,5 @@
 program example_FE_2D
-   use linear_pack, only: general_square_matrix
-   use quadrature_module, only: wp => quadrature_wp
+   use mod_kinds
    use mod_FE_2D
 
    implicit none
@@ -26,11 +25,6 @@ program example_FE_2D
    real(wp), allocatable :: M(:, :)
    integer, allocatable :: T(:, :)
 
-   real(wp), allocatable :: solution(:, :)
-
-   integer :: N_basis
-   real(wp) :: h_basis, max_FE_error
-
    integer :: index
 
    print *, "------------------start------------------"

example/example_2D/example_field_2D.f90

@@ -1,6 +1,5 @@
 program example_field_2D
-   use linear_pack, only: general_square_matrix
-   use quadrature_module, only: wp => quadrature_wp
+   use mod_kinds
    use mod_FE_2D
 
    implicit none
@@ -10,26 +9,26 @@ program example_field_2D
 
    !> Field parameter input.
    real(wp) :: left = 0.0_wp, right = 1.0_wp, bottom = 0.0_wp, top = 1.0_wp
-      !! The problem domain is [left,right]*[bottom,top].
+   !! The problem domain is [left,right]*[bottom,top].
    integer :: Nh_parition = 2, Nv_parition = 2
-      !! The number of partition in horizontal and vertical direction.
+   !! The number of partition in horizontal and vertical direction.
    integer :: Nh_basis, Nv_basis
-      !! The number of FE basis functions in horizontal and vertical direction.
+   !! The number of FE basis functions in horizontal and vertical direction.
    integer :: Gauss_point_number = 8
-      !! The number of Gauss Quadrature points.
+   !! The number of Gauss Quadrature points.
    !> basis_type: the type of the FE.
    integer :: basis_type = 201
-      !! 201: 2D linear
-      !! 202: 2D quadratic
+   !! 201: 2D linear
+   !! 202: 2D quadratic
    integer :: mesh_type = 503
-      !! 503: 2D triangular mesh;
-      !! 504: 2D quadrilateral mesh;
+   !! 503: 2D triangular mesh;
+   !! 504: 2D quadrilateral mesh;
 
    select case (basis_type)
-   case (201)
+    case (201)
       Nh_basis = Nh_parition
       Nv_basis = Nv_parition
-   case (202)
+    case (202)
       Nh_basis = 2*Nh_parition
       Nv_basis = 2*Nv_parition
    end select
@@ -38,13 +37,13 @@ program example_field_2D
    !! TODO: Input the field information.
    !> Initialize and check the field information.
    call field_info%init(left, right, bottom, top, &
-                        Nh_parition, Nv_parition, &
-                        Nh_basis, Nv_basis, &
-                        Gauss_point_number=Gauss_point_number, &
-                        tolerance=1.0e-6_wp, &
-                        trial_basis_type=basis_type, &
-                        test_basis_type=basis_type, &
-                        mesh_type=mesh_type)
+      Nh_parition, Nv_parition, &
+      Nh_basis, Nv_basis, &
+      Gauss_point_number=Gauss_point_number, &
+      tolerance=1.0e-6_wp, &
+      trial_basis_type=basis_type, &
+      test_basis_type=basis_type, &
+      mesh_type=mesh_type)
 
    !> Use the format (A,2F6.3) to print the field information.
    call field_info%print('(A,2F6.3)')

example/example_2D/example_generate_2D.f90

@@ -1,29 +1,28 @@
 program example_generate_2D
-   use linear_pack, only: general_square_matrix
-   use quadrature_module, only: wp => quadrature_wp
+   use mod_kinds
    use mod_FE_2D
-
+   
    implicit none
 
    !> input information of field
    type(field) :: field_info
 
    !> Field parameter input.
    real(wp) :: left = 0.0_wp, right = 1.0_wp, bottom = 0.0_wp, top = 1.0_wp
-      !! The problem domain is [left,right]*[bottom,top].
+   !! The problem domain is [left,right]*[bottom,top].
    integer :: Nh_parition = 2, Nv_parition = 2
-      !! The number of partition in horizontal and vertical direction.
+   !! The number of partition in horizontal and vertical direction.
    integer :: Nh_basis, Nv_basis
-      !! The number of FE basis functions in horizontal and vertical direction.
+   !! The number of FE basis functions in horizontal and vertical direction.
    integer :: Gauss_point_number = 8
-      !! The number of Gauss Quadrature points.
+   !! The number of Gauss Quadrature points.
    !> basis_type: the type of the FE.
    integer :: basis_type = 202
-      !! 201: 2D linear
-      !! 202: 2D quadratic
+   !! 201: 2D linear
+   !! 202: 2D quadratic
    integer :: mesh_type = 504
-      !! 503: 2D triangular mesh;
-      !! 504: 2D quadrilateral mesh;
+   !! 503: 2D triangular mesh;
+   !! 504: 2D quadrilateral mesh;
 
    real(wp), allocatable :: M(:, :)
    integer, allocatable :: T(:, :)
@@ -33,24 +32,24 @@ program example_generate_2D
    integer :: index, index_2
 
    select case (basis_type)
-   case (201)
+    case (201)
       Nh_basis = Nh_parition
       Nv_basis = Nv_parition
-   case (202)
+    case (202)
       Nh_basis = 2*Nh_parition
       Nv_basis = 2*Nv_parition
    end select
 
    print *, "------------------Information Matrix Generate------------------"
 
    call field_info%init(left, right, bottom, top, &
-                        Nh_parition, Nv_parition, &
-                        Nh_basis, Nv_basis, &
-                        Gauss_point_number=Gauss_point_number, &
-                        tolerance=1.0e-6_wp, &
-                        trial_basis_type=basis_type, &
-                        test_basis_type=basis_type, &
-                        mesh_type=mesh_type)
+      Nh_parition, Nv_parition, &
+      Nh_basis, Nv_basis, &
+      Gauss_point_number=Gauss_point_number, &
+      tolerance=1.0e-6_wp, &
+      trial_basis_type=basis_type, &
+      test_basis_type=basis_type, &
+      mesh_type=mesh_type)
    !> Check the field information.
    call field_info%check()
    !> Print the field information.

fpm.toml

@@ -20,5 +20,7 @@ source-form = "free"
 M_attr = { git = "https://github.com/urbanjost/M_attr.git" }
 # M_CLI2 = { git = "https://github.com/urbanjost/M_CLI2.git" }
 # toml-f = { git = "https://github.com/toml-f/toml-f.git" }
+LSQR = { git = "https://github.com/jacobwilliams/lsqr.git" }
 # test-drive = { git = "https://github.com/fortran-lang/test-drive.git" }
 # fortran-lua53 = { git = "https://github.com/interkosmos/fortran-lua53.git" }
+