poisson_iga_solver_3d Program

Uses

  • program~~poisson_iga_solver_3d~~UsesGraph program~poisson_iga_solver_3d poisson_iga_solver_3d fortime fortime program~poisson_iga_solver_3d->fortime module~forcad forcad program~poisson_iga_solver_3d->module~forcad module~forcad_utils forcad_utils program~poisson_iga_solver_3d->module~forcad_utils module~forcad_kinds forcad_kinds module~forcad->module~forcad_kinds module~forcad_nurbs_curve forcad_nurbs_curve module~forcad->module~forcad_nurbs_curve module~forcad_nurbs_surface forcad_nurbs_surface module~forcad->module~forcad_nurbs_surface module~forcad_nurbs_volume forcad_nurbs_volume module~forcad->module~forcad_nurbs_volume module~forcad_utils->module~forcad_kinds module~forcad_nurbs_curve->module~forcad_utils module~forcad_nurbs_curve->module~forcad_kinds module~forcad_nurbs_surface->module~forcad_utils module~forcad_nurbs_surface->module~forcad_kinds module~forcad_nurbs_volume->module~forcad_utils module~forcad_nurbs_volume->module~forcad_kinds

Solves the 3D Poisson problem using Isogeometric Analysis (IGA).

This code solves the equation: using a B-spline volume on a rectangular domain.

The solution is discretized using tensor-product B-spline basis functions over a structured control net. The global linear system is assembled and solved using a basic internal Cholesky solver.

The resulting solution is exported to VTK format for visualization.

The L2 error norm with respect to the exact solution is computed and printed.

Note

This implementation uses B-spline geometry (no rational weights), hence the volume is not a full NURBS.

Warning

"Slow solver" The solver uses an internal Cholesky factorization which is not optimized. For large systems, consider replacing it with a more scalable external solver.

Domain size and number of control points Number knots to insert in each direction Mode numbers for the source term and exact solution Resolution of the visualization grid filename for VTK export

Construct the NURBS surface For simplicity, set_hexahedron creates a rectangular surface with uniform knot spacing For more complex geometries, use surf%set() with knots, continuity,...

Insert knots in the first and second directions Extract geometry and mesh structure Assemble global stiffness matrix and load vector Apply homogeneous Dirichlet boundary conditions Solve the linear system K·X = b Export solution at control points to VTK

Interpolate solution and export field Compute the L2 error norm


Calls

program~~poisson_iga_solver_3d~~CallsGraph program~poisson_iga_solver_3d poisson_iga_solver_3d none~basis~3 nurbs_volume%basis program~poisson_iga_solver_3d->none~basis~3 none~get_degree~3 nurbs_volume%get_degree program~poisson_iga_solver_3d->none~get_degree~3 none~get_nc~3 nurbs_volume%get_nc program~poisson_iga_solver_3d->none~get_nc~3 none~get_xc~3 nurbs_volume%get_Xc program~poisson_iga_solver_3d->none~get_xc~3 proc~ansatz~3 nurbs_volume%ansatz program~poisson_iga_solver_3d->proc~ansatz~3 proc~cmp_elem~3 nurbs_volume%cmp_elem program~poisson_iga_solver_3d->proc~cmp_elem~3 proc~create~3 nurbs_volume%create program~poisson_iga_solver_3d->proc~create~3 proc~exact_solution exact_solution program~poisson_iga_solver_3d->proc~exact_solution proc~export_xc~3 nurbs_volume%export_Xc program~poisson_iga_solver_3d->proc~export_xc~3 proc~export_xg~3 nurbs_volume%export_Xg program~poisson_iga_solver_3d->proc~export_xg~3 proc~finalize~3 nurbs_volume%finalize program~poisson_iga_solver_3d->proc~finalize~3 proc~insert_knots~3 nurbs_volume%insert_knots program~poisson_iga_solver_3d->proc~insert_knots~3 proc~set_hexahedron nurbs_volume%set_hexahedron program~poisson_iga_solver_3d->proc~set_hexahedron proc~solve solve program~poisson_iga_solver_3d->proc~solve proc~source_term source_term program~poisson_iga_solver_3d->proc~source_term timer_start timer_start program~poisson_iga_solver_3d->timer_start timer_stop timer_stop program~poisson_iga_solver_3d->timer_stop proc~basis_scalar~3 nurbs_volume%basis_scalar none~basis~3->proc~basis_scalar~3 proc~basis_vector~3 nurbs_volume%basis_vector none~basis~3->proc~basis_vector~3 proc~get_degree_all~2 nurbs_volume%get_degree_all none~get_degree~3->proc~get_degree_all~2 proc~get_degree_dir~2 nurbs_volume%get_degree_dir none~get_degree~3->proc~get_degree_dir~2 proc~get_nc_all~2 nurbs_volume%get_nc_all none~get_nc~3->proc~get_nc_all~2 proc~get_nc_dir~2 nurbs_volume%get_nc_dir none~get_nc~3->proc~get_nc_dir~2 proc~get_xc_all~3 nurbs_volume%get_Xc_all none~get_xc~3->proc~get_xc_all~3 proc~get_xcid~3 nurbs_volume%get_Xcid none~get_xc~3->proc~get_xcid~3 proc~get_xci~3 nurbs_volume%get_Xci none~get_xc~3->proc~get_xci~3 proc~ansatz~3->proc~cmp_elem~3 interface~gauss_leg gauss_leg proc~ansatz~3->interface~gauss_leg interface~ndgrid ndgrid proc~ansatz~3->interface~ndgrid interface~unique unique proc~ansatz~3->interface~unique none~derivative~3 nurbs_volume%derivative proc~ansatz~3->none~derivative~3 none~set~3 nurbs_volume%set proc~ansatz~3->none~set~3 proc~det det proc~ansatz~3->proc~det proc~inv inv proc~ansatz~3->proc~inv interface~elemconn_cn elemConn_Cn proc~cmp_elem~3->interface~elemconn_cn proc~cmp_elem~3->interface~unique proc~get_multiplicity~3 nurbs_volume%get_multiplicity proc~cmp_elem~3->proc~get_multiplicity~3 interface~compute_xg~3 compute_Xg proc~create~3->interface~compute_xg~3 proc~create~3->interface~ndgrid proc~is_rational~3 nurbs_volume%is_rational proc~create~3->proc~is_rational~3 proc~cmp_elem_xc_vis~3 nurbs_volume%cmp_elem_Xc_vis proc~export_xc~3->proc~cmp_elem_xc_vis~3 proc~export_vtk_legacy export_vtk_legacy proc~export_xc~3->proc~export_vtk_legacy proc~cmp_elem_xg_vis~3 nurbs_volume%cmp_elem_Xg_vis proc~export_xg~3->proc~cmp_elem_xg_vis~3 proc~export_xg~3->proc~export_vtk_legacy interface~compute_multiplicity compute_multiplicity proc~insert_knots~3->interface~compute_multiplicity none~get_knot~3 nurbs_volume%get_knot proc~insert_knots~3->none~get_knot~3 proc~insert_knots~3->none~set~3 proc~findspan findspan proc~insert_knots~3->proc~findspan proc~insert_knot_a_5_1 insert_knot_A_5_1 proc~insert_knots~3->proc~insert_knot_a_5_1 proc~set_hexahedron->none~set~3 proc~hexahedron_xc hexahedron_Xc proc~set_hexahedron->proc~hexahedron_xc proc~compute_multiplicity1 compute_multiplicity1 interface~compute_multiplicity->proc~compute_multiplicity1 proc~compute_multiplicity2 compute_multiplicity2 interface~compute_multiplicity->proc~compute_multiplicity2 proc~compute_xg_bspline_3d compute_Xg_bspline_3d interface~compute_xg~3->proc~compute_xg_bspline_3d proc~compute_xg_bspline_3d_1point compute_Xg_bspline_3d_1point interface~compute_xg~3->proc~compute_xg_bspline_3d_1point proc~compute_xg_nurbs_3d compute_Xg_nurbs_3d interface~compute_xg~3->proc~compute_xg_nurbs_3d proc~compute_xg_nurbs_3d_1point compute_Xg_nurbs_3d_1point interface~compute_xg~3->proc~compute_xg_nurbs_3d_1point proc~cmp_elemconn_cn_l cmp_elemConn_Cn_L interface~elemconn_cn->proc~cmp_elemconn_cn_l proc~cmp_elemconn_cn_s cmp_elemConn_Cn_S interface~elemconn_cn->proc~cmp_elemconn_cn_s proc~cmp_elemconn_cn_v cmp_elemConn_Cn_V interface~elemconn_cn->proc~cmp_elemconn_cn_v proc~gauss_legendre_1d gauss_legendre_1D interface~gauss_leg->proc~gauss_legendre_1d proc~gauss_legendre_2d gauss_legendre_2D interface~gauss_leg->proc~gauss_legendre_2d proc~gauss_legendre_3d gauss_legendre_3D interface~gauss_leg->proc~gauss_legendre_3d proc~ndgrid2 ndgrid2 interface~ndgrid->proc~ndgrid2 proc~ndgrid3 ndgrid3 interface~ndgrid->proc~ndgrid3 proc~unique_integer unique_integer interface~unique->proc~unique_integer proc~unique_real unique_real interface~unique->proc~unique_real proc~derivative_scalar~3 nurbs_volume%derivative_scalar none~derivative~3->proc~derivative_scalar~3 proc~derivative_vector~3 nurbs_volume%derivative_vector none~derivative~3->proc~derivative_vector~3 proc~get_knot_all~3 nurbs_volume%get_knot_all none~get_knot~3->proc~get_knot_all~3 proc~get_knoti~3 nurbs_volume%get_knoti none~get_knot~3->proc~get_knoti~3 proc~set1~3 nurbs_volume%set1 none~set~3->proc~set1~3 proc~set2~3 nurbs_volume%set2 none~set~3->proc~set2~3 proc~set3~3 nurbs_volume%set3 none~set~3->proc~set3~3 proc~set4~3 nurbs_volume%set4 none~set~3->proc~set4~3 proc~basis_scalar~3->proc~is_rational~3 interface~compute_tgc~3 compute_Tgc proc~basis_scalar~3->interface~compute_tgc~3 proc~basis_vector~3->interface~ndgrid proc~basis_vector~3->proc~is_rational~3 proc~basis_vector~3->interface~compute_tgc~3 interface~elemconn_c0 elemConn_C0 proc~cmp_elem_xc_vis~3->interface~elemconn_c0 proc~cmp_elem_xg_vis~3->interface~elemconn_c0 proc~get_multiplicity~3->interface~compute_multiplicity proc~get_nc_dir~2->interface~compute_multiplicity proc~inv->proc~det proc~inv->proc~inv proc~compute_tgc_bspline_3d_scalar compute_Tgc_bspline_3d_scalar interface~compute_tgc~3->proc~compute_tgc_bspline_3d_scalar proc~compute_tgc_bspline_3d_vector compute_Tgc_bspline_3d_vector interface~compute_tgc~3->proc~compute_tgc_bspline_3d_vector proc~compute_tgc_nurbs_3d_scalar compute_Tgc_nurbs_3d_scalar interface~compute_tgc~3->proc~compute_tgc_nurbs_3d_scalar proc~compute_tgc_nurbs_3d_vector compute_Tgc_nurbs_3d_vector interface~compute_tgc~3->proc~compute_tgc_nurbs_3d_vector proc~cmp_elemconn_c0_l cmp_elemConn_C0_L interface~elemconn_c0->proc~cmp_elemconn_c0_l proc~cmp_elemconn_c0_s cmp_elemConn_C0_S interface~elemconn_c0->proc~cmp_elemconn_c0_s proc~cmp_elemconn_c0_v cmp_elemConn_C0_V interface~elemconn_c0->proc~cmp_elemconn_c0_v proc~basis_bspline basis_bspline proc~compute_xg_bspline_3d->proc~basis_bspline proc~kron kron proc~compute_xg_bspline_3d->proc~kron proc~compute_xg_bspline_3d_1point->proc~basis_bspline proc~compute_xg_bspline_3d_1point->proc~kron proc~cmp_tgc_3d cmp_Tgc_3d proc~compute_xg_nurbs_3d->proc~cmp_tgc_3d proc~compute_xg_nurbs_3d_1point->proc~basis_bspline proc~compute_xg_nurbs_3d_1point->proc~kron proc~derivative_scalar~3->proc~is_rational~3 interface~compute_dtgc~3 compute_dTgc proc~derivative_scalar~3->interface~compute_dtgc~3 proc~derivative_vector~3->interface~ndgrid proc~derivative_vector~3->proc~is_rational~3 proc~derivative_vector~3->interface~compute_dtgc~3 proc~gauss_legendre gauss_legendre proc~gauss_legendre_1d->proc~gauss_legendre proc~gauss_legendre_2d->interface~ndgrid proc~gauss_legendre_2d->proc~gauss_legendre proc~gauss_legendre_2d->proc~kron proc~gauss_legendre_3d->interface~ndgrid proc~gauss_legendre_3d->proc~gauss_legendre proc~gauss_legendre_3d->proc~kron proc~cmp_degree~3 nurbs_volume%cmp_degree proc~set1~3->proc~cmp_degree~3 proc~cmp_nc~3 nurbs_volume%cmp_nc proc~set1~3->proc~cmp_nc~3 proc~set2~3->proc~cmp_nc~3 proc~compute_knot_vector compute_knot_vector proc~set2~3->proc~compute_knot_vector proc~set3~3->proc~cmp_degree~3 proc~compute_dtgc_bspline_3d_scalar compute_dTgc_bspline_3d_scalar interface~compute_dtgc~3->proc~compute_dtgc_bspline_3d_scalar proc~compute_dtgc_bspline_3d_vector compute_dTgc_bspline_3d_vector interface~compute_dtgc~3->proc~compute_dtgc_bspline_3d_vector proc~compute_dtgc_nurbs_3d_scalar compute_dTgc_nurbs_3d_scalar interface~compute_dtgc~3->proc~compute_dtgc_nurbs_3d_scalar proc~compute_dtgc_nurbs_3d_vector compute_dTgc_nurbs_3d_vector interface~compute_dtgc~3->proc~compute_dtgc_nurbs_3d_vector proc~cmp_degree~3->proc~get_multiplicity~3 proc~cmp_nc~3->interface~compute_multiplicity proc~cmp_tgc_3d->proc~basis_bspline proc~cmp_tgc_3d->proc~kron proc~repelem repelem proc~compute_knot_vector->proc~repelem proc~compute_tgc_bspline_3d_scalar->proc~basis_bspline proc~compute_tgc_bspline_3d_scalar->proc~kron proc~compute_tgc_bspline_3d_vector->proc~basis_bspline proc~compute_tgc_bspline_3d_vector->proc~kron proc~compute_tgc_nurbs_3d_scalar->proc~basis_bspline proc~compute_tgc_nurbs_3d_scalar->proc~kron proc~compute_tgc_nurbs_3d_vector->proc~basis_bspline proc~compute_tgc_nurbs_3d_vector->proc~kron proc~compute_dtgc_bspline_3d_scalar->proc~kron interface~basis_bspline_der basis_bspline_der proc~compute_dtgc_bspline_3d_scalar->interface~basis_bspline_der proc~compute_dtgc_bspline_3d_vector->proc~kron proc~compute_dtgc_bspline_3d_vector->interface~basis_bspline_der proc~compute_dtgc_nurbs_3d_scalar->proc~kron proc~compute_dtgc_nurbs_3d_scalar->interface~basis_bspline_der proc~compute_dtgc_nurbs_3d_vector->proc~kron proc~compute_dtgc_nurbs_3d_vector->interface~basis_bspline_der proc~basis_bspline_der_a basis_bspline_der_A interface~basis_bspline_der->proc~basis_bspline_der_a proc~basis_bspline_der_b basis_bspline_der_B interface~basis_bspline_der->proc~basis_bspline_der_b

Variables

Type Attributes Name Initial
real(kind=rk), allocatable :: K(:,:)

Global stiffness matrix

real(kind=rk) :: L(3)

Domain size

real(kind=rk), allocatable :: T(:)

Basis function values at quadrature point
real(kind=rk), allocatable :: X(:,:)

Global solution vector

real(kind=rk), allocatable :: Xc(:,:)

Control point coordinates
real(kind=rk), allocatable :: Xg(:)

Physical coordinates at quadrature point
real(kind=rk), allocatable :: b(:)

Global right-hand side vector

real(kind=rk), allocatable :: dT(:,:)

Derivatives of basis functions at quadrature point
real(kind=rk) :: dV

Differential element volume

integer, allocatable :: dirichlet_id(:)

Indices for Dirichlet boundary conditions
integer :: dof

Degrees of freedom per control point (1 for scalar field)
integer, allocatable :: elem(:,:)

Element connectivity matrix
integer, allocatable :: elem_e(:)

Local connectivity of current element
character(len=256) :: filename

Filename for VTK export
integer :: i

Generic loop index
integer :: ie

Element index
integer :: ig

Quadrature (Gauss) point index
integer :: ki(3)

Number of knots to insert in each direction
real(kind=rk) :: l2_error

L2 error norm accumulator

integer :: m(3)

Mode numbers for source and exact solution
integer :: nc(3)

Number of control points in each direction
integer :: nct

Total number of control points
integer :: ndof

Total degrees of freedom
integer :: nelem

Number of elements
integer :: nnelem

Number of local nodes per element
real(kind=rk), parameter :: pi = acos(-1.0_rk)

Constant

integer :: res(3)

Visualization resolution
type(nurbs_volume) :: surf

NURBS volume object
type(timer) :: ti

Timer object for performance measurement
real(kind=rk), allocatable :: u_h(:,:)

Interpolated solution field on grid

Functions

pure function exact_solution(p, d, n) result(u)

Computes the exact solution corresponding to the source term

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: p(3)

Coordinates (x, y, z)
real(kind=rk), intent(in) :: d(3)

Domain size (L1, L2, L3)
integer, intent(in) :: n(3)

Mode numbers (m1, m2, m3)

Return Value real(kind=rk)

pure function source_term(p, d, n) result(f)

Computes the source function

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: p(3)

Coordinates (x, y, z)
real(kind=rk), intent(in) :: d(3)

Domain size (L1, L2, L3)
integer, intent(in) :: n(3)

Mode numbers (m1, m2, m3)

Return Value real(kind=rk)