construct_square_first Module Function

    module function construct_square_first(id, global_coordinate, cell_info, integration) result(element)
        implicit none
        integer(int32), intent(in) :: id
        type(type_dp_3d), pointer, intent(in) :: global_coordinate
        type(type_vtk_cell), intent(in) :: cell_info
        type(type_geometry_settings), intent(in) :: integration
        class(abst_element), allocatable :: element

        integer(int32) :: i

        if (allocated(element)) deallocate (element)
        allocate (type_square_second :: element)

        element%id = id
        element%type = cell_info%cell_type
        element%group = cell_info%cell_entity_id
        element%dimension = cell_info%get_dimension()
        element%order = cell_info%get_order()

        element%num_nodes = cell_info%num_nodes_in_cell
        allocate (element%connectivity(element%num_nodes))
        element%connectivity(:) = cell_info%connectivity(1:element%num_nodes)

        allocate (element%x(element%num_nodes))
        allocate (element%y(element%num_nodes))
        allocate (element%z(element%num_nodes))
        do i = 1, element%num_nodes
            nullify (element%x(i)%val)
            nullify (element%y(i)%val)
            nullify (element%z(i)%val)
            element%x(i)%val => global_coordinate%x(element%connectivity(i))
            element%y(i)%val => global_coordinate%y(element%connectivity(i))
            element%z(i)%val => global_coordinate%z(element%connectivity(i))
        end do

        select case (integration%integration_type)
        case ("full")
            element%num_gauss = 4
            call allocate_array(element%weight, element%num_gauss)
            call allocate_array(element%gauss, element%dimension, element%num_gauss)

            element%weight(:) = [1.0d0, 1.0d0, 1.0d0, 1.0d0]
            element%gauss(:, 1) = [-sqrt(1.0d0 / 3.0d0), -sqrt(1.0d0 / 3.0d0)]
            element%gauss(:, 2) = [-sqrt(1.0d0 / 3.0d0), sqrt(1.0d0 / 3.0d0)]
            element%gauss(:, 3) = [sqrt(1.0d0 / 3.0d0), sqrt(1.0d0 / 3.0d0)]
            element%gauss(:, 4) = [sqrt(1.0d0 / 3.0d0), -sqrt(1.0d0 / 3.0d0)]
        case ("reduced")
            element%num_gauss = 1_int32
            call allocate_array(element%weight, element%num_gauss)
            call allocate_array(element%gauss, element%dimension, element%num_gauss)

            element%weight(:) = [4.0d0]
            element%gauss(:, 1) = [0.0d0, 0.0d0]
        case ("free")
            element%num_gauss = 4
            call allocate_array(element%weight, element%num_gauss)
            call allocate_array(element%gauss, element%dimension, element%num_gauss)

            element%weight(:) = [1.0d0, 1.0d0, 1.0d0, 1.0d0]
            element%gauss(:, 1) = [-integration%integration_points, -integration%integration_points]
            element%gauss(:, 2) = [-integration%integration_points, integration%integration_points]
            element%gauss(:, 3) = [integration%integration_points, integration%integration_points]
            element%gauss(:, 4) = [integration%integration_points, -integration%integration_points]
        end select

        if (associated(element%interpolate)) nullify (element%interpolate)
        element%interpolate => interpolate

        if (associated(element%get_connectivity)) nullify (element%get_connectivity)
        element%get_connectivity => get_connectivity

    end function construct_square_first