Hello everyone!
I'm working on a 1D Finite Element Method (FEM) solver that supports p-version formulations. The p-version uses hierarchical shape functions based on Legendre polynomials to increase the polynomial order of approximation within each element. The problem involves a bar fixed at both ends and subjected to a smooth distributed body force. While implementing the p-FEM approach, I observed that although the Dirichlet boundary conditions are correctly enforced (the displacements at the endpoints are zero), one of the internal displacement values becomes unexpectedly negative, which does not make physical sense given the nature of the body force. The stiffness matrix appears correctly assembled, symmetric, and positive-definite, and the body force integration follows the expected Gaussian quadrature process. I suspect the issue may lie in how the fictitious nodal coordinates are constructed within each element in the p-version, possibly leading to incorrect evaluation of the body force or shape function derivatives. Despite checking the Jacobian signs and shape function values, the resulting solution still shows a large negative value internally, which suggests a deeper inconsistency in either the mapping or the integration process.
Can somebody help me?