poisson_iga_solver_2d Program

Uses

  • program~~poisson_iga_solver_2d~~UsesGraph program~poisson_iga_solver_2d poisson_iga_solver_2d fortime fortime program~poisson_iga_solver_2d->fortime module~forcad forcad program~poisson_iga_solver_2d->module~forcad module~forcad_utils forcad_utils program~poisson_iga_solver_2d->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 2D Poisson problem using Isogeometric Analysis (IGA).

This code solves the equation: using a B-spline surface 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 surface 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_tetragon 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_2d~~CallsGraph program~poisson_iga_solver_2d poisson_iga_solver_2d none~basis nurbs_surface%basis program~poisson_iga_solver_2d->none~basis none~get_degree nurbs_surface%get_degree program~poisson_iga_solver_2d->none~get_degree none~get_nc nurbs_surface%get_nc program~poisson_iga_solver_2d->none~get_nc none~get_xc nurbs_surface%get_Xc program~poisson_iga_solver_2d->none~get_xc proc~ansatz nurbs_surface%ansatz program~poisson_iga_solver_2d->proc~ansatz proc~cmp_elem nurbs_surface%cmp_elem program~poisson_iga_solver_2d->proc~cmp_elem proc~create nurbs_surface%create program~poisson_iga_solver_2d->proc~create proc~exact_solution~2 exact_solution program~poisson_iga_solver_2d->proc~exact_solution~2 proc~export_xc nurbs_surface%export_Xc program~poisson_iga_solver_2d->proc~export_xc proc~export_xg nurbs_surface%export_Xg program~poisson_iga_solver_2d->proc~export_xg proc~finalize nurbs_surface%finalize program~poisson_iga_solver_2d->proc~finalize proc~insert_knots nurbs_surface%insert_knots program~poisson_iga_solver_2d->proc~insert_knots proc~set_tetragon nurbs_surface%set_tetragon program~poisson_iga_solver_2d->proc~set_tetragon proc~solve solve program~poisson_iga_solver_2d->proc~solve proc~source_term~2 source_term program~poisson_iga_solver_2d->proc~source_term~2 timer_start timer_start program~poisson_iga_solver_2d->timer_start timer_stop timer_stop program~poisson_iga_solver_2d->timer_stop proc~basis_scalar nurbs_surface%basis_scalar none~basis->proc~basis_scalar proc~basis_vector nurbs_surface%basis_vector none~basis->proc~basis_vector proc~get_degree_all nurbs_surface%get_degree_all none~get_degree->proc~get_degree_all proc~get_degree_dir nurbs_surface%get_degree_dir none~get_degree->proc~get_degree_dir proc~get_nc_all nurbs_surface%get_nc_all none~get_nc->proc~get_nc_all proc~get_nc_dir nurbs_surface%get_nc_dir none~get_nc->proc~get_nc_dir proc~get_xc_all nurbs_surface%get_Xc_all none~get_xc->proc~get_xc_all proc~get_xci nurbs_surface%get_Xci none~get_xc->proc~get_xci proc~get_xcid nurbs_surface%get_Xcid none~get_xc->proc~get_xcid proc~ansatz->proc~cmp_elem interface~gauss_leg gauss_leg proc~ansatz->interface~gauss_leg interface~ndgrid ndgrid proc~ansatz->interface~ndgrid interface~unique unique proc~ansatz->interface~unique none~derivative nurbs_surface%derivative proc~ansatz->none~derivative none~set nurbs_surface%set proc~ansatz->none~set proc~det det proc~ansatz->proc~det proc~inv inv proc~ansatz->proc~inv interface~elemconn_cn elemConn_Cn proc~cmp_elem->interface~elemconn_cn proc~cmp_elem->interface~unique proc~get_multiplicity nurbs_surface%get_multiplicity proc~cmp_elem->proc~get_multiplicity interface~compute_xg compute_Xg proc~create->interface~compute_xg proc~create->interface~ndgrid proc~is_rational nurbs_surface%is_rational proc~create->proc~is_rational proc~cmp_elem_xc_vis nurbs_surface%cmp_elem_Xc_vis proc~export_xc->proc~cmp_elem_xc_vis proc~export_vtk_legacy export_vtk_legacy proc~export_xc->proc~export_vtk_legacy proc~cmp_elem_xg_vis nurbs_surface%cmp_elem_Xg_vis proc~export_xg->proc~cmp_elem_xg_vis proc~export_xg->proc~export_vtk_legacy interface~compute_multiplicity compute_multiplicity proc~insert_knots->interface~compute_multiplicity none~get_knot nurbs_surface%get_knot proc~insert_knots->none~get_knot proc~insert_knots->none~set proc~findspan findspan proc~insert_knots->proc~findspan proc~insert_knot_a_5_1 insert_knot_A_5_1 proc~insert_knots->proc~insert_knot_a_5_1 proc~insert_knots->proc~is_rational proc~set_tetragon->none~set proc~tetragon_xc tetragon_Xc proc~set_tetragon->proc~tetragon_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_2d compute_Xg_bspline_2d interface~compute_xg->proc~compute_xg_bspline_2d proc~compute_xg_bspline_2d_1point compute_Xg_bspline_2d_1point interface~compute_xg->proc~compute_xg_bspline_2d_1point proc~compute_xg_nurbs_2d compute_Xg_nurbs_2d interface~compute_xg->proc~compute_xg_nurbs_2d proc~compute_xg_nurbs_2d_1point compute_Xg_nurbs_2d_1point interface~compute_xg->proc~compute_xg_nurbs_2d_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 nurbs_surface%derivative_scalar none~derivative->proc~derivative_scalar proc~derivative_vector nurbs_surface%derivative_vector none~derivative->proc~derivative_vector proc~get_knot_all nurbs_surface%get_knot_all none~get_knot->proc~get_knot_all proc~get_knoti nurbs_surface%get_knoti none~get_knot->proc~get_knoti proc~set1 nurbs_surface%set1 none~set->proc~set1 proc~set2 nurbs_surface%set2 none~set->proc~set2 proc~set3 nurbs_surface%set3 none~set->proc~set3 proc~set4 nurbs_surface%set4 none~set->proc~set4 proc~basis_scalar->proc~is_rational interface~compute_tgc compute_Tgc proc~basis_scalar->interface~compute_tgc proc~basis_vector->interface~ndgrid proc~basis_vector->proc~is_rational proc~basis_vector->interface~compute_tgc interface~elemconn_c0 elemConn_C0 proc~cmp_elem_xc_vis->interface~elemconn_c0 proc~cmp_elem_xg_vis->interface~elemconn_c0 proc~get_multiplicity->interface~compute_multiplicity proc~get_nc_dir->interface~compute_multiplicity proc~inv->proc~det proc~inv->proc~inv proc~compute_tgc_bspline_2d_scalar compute_Tgc_bspline_2d_scalar interface~compute_tgc->proc~compute_tgc_bspline_2d_scalar proc~compute_tgc_bspline_2d_vector compute_Tgc_bspline_2d_vector interface~compute_tgc->proc~compute_tgc_bspline_2d_vector proc~compute_tgc_nurbs_2d_scalar compute_Tgc_nurbs_2d_scalar interface~compute_tgc->proc~compute_tgc_nurbs_2d_scalar proc~compute_tgc_nurbs_2d_vector compute_Tgc_nurbs_2d_vector interface~compute_tgc->proc~compute_tgc_nurbs_2d_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_2d->proc~basis_bspline proc~kron kron proc~compute_xg_bspline_2d->proc~kron proc~compute_xg_bspline_2d_1point->proc~basis_bspline proc~compute_xg_bspline_2d_1point->proc~kron proc~cmp_tgc_2d cmp_Tgc_2d proc~compute_xg_nurbs_2d->proc~cmp_tgc_2d proc~compute_xg_nurbs_2d_1point->proc~basis_bspline proc~compute_xg_nurbs_2d_1point->proc~kron proc~derivative_scalar->proc~is_rational interface~compute_dtgc compute_dTgc proc~derivative_scalar->interface~compute_dtgc proc~derivative_vector->interface~ndgrid proc~derivative_vector->proc~is_rational proc~derivative_vector->interface~compute_dtgc 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 nurbs_surface%cmp_degree proc~set1->proc~cmp_degree proc~cmp_nc nurbs_surface%cmp_nc proc~set1->proc~cmp_nc proc~set2->proc~cmp_nc proc~compute_knot_vector compute_knot_vector proc~set2->proc~compute_knot_vector proc~set3->proc~cmp_degree proc~compute_dtgc_bspline_2d_scalar compute_dTgc_bspline_2d_scalar interface~compute_dtgc->proc~compute_dtgc_bspline_2d_scalar proc~compute_dtgc_bspline_2d_vector compute_dTgc_bspline_2d_vector interface~compute_dtgc->proc~compute_dtgc_bspline_2d_vector proc~compute_dtgc_nurbs_2d_scalar compute_dTgc_nurbs_2d_scalar interface~compute_dtgc->proc~compute_dtgc_nurbs_2d_scalar proc~compute_dtgc_nurbs_2d_vector compute_dTgc_nurbs_2d_vector interface~compute_dtgc->proc~compute_dtgc_nurbs_2d_vector proc~cmp_degree->proc~get_multiplicity proc~cmp_nc->interface~compute_multiplicity proc~cmp_tgc_2d->proc~basis_bspline proc~cmp_tgc_2d->proc~kron proc~repelem repelem proc~compute_knot_vector->proc~repelem proc~compute_tgc_bspline_2d_scalar->proc~basis_bspline proc~compute_tgc_bspline_2d_scalar->proc~kron proc~compute_tgc_bspline_2d_vector->proc~basis_bspline proc~compute_tgc_bspline_2d_vector->proc~kron proc~compute_tgc_nurbs_2d_scalar->proc~basis_bspline proc~compute_tgc_nurbs_2d_scalar->proc~kron proc~compute_tgc_nurbs_2d_vector->proc~basis_bspline proc~compute_tgc_nurbs_2d_vector->proc~kron proc~compute_dtgc_bspline_2d_scalar->proc~kron interface~basis_bspline_der basis_bspline_der proc~compute_dtgc_bspline_2d_scalar->interface~basis_bspline_der proc~compute_dtgc_bspline_2d_vector->proc~kron proc~compute_dtgc_bspline_2d_vector->interface~basis_bspline_der proc~compute_dtgc_nurbs_2d_scalar->proc~kron proc~compute_dtgc_nurbs_2d_scalar->interface~basis_bspline_der proc~compute_dtgc_nurbs_2d_vector->proc~kron proc~compute_dtgc_nurbs_2d_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(2)

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) :: dA

Differential element area

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

Derivatives of basis functions at quadrature point
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(2)

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

L2 error norm accumulator

integer :: m(2)

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

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(2)

Visualization resolution
type(nurbs_surface) :: surf

NURBS surface 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(2)

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

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

Mode numbers (m1, m2)

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(2)

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

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

Mode numbers (m1, m2)

Return Value real(kind=rk)