I am mechanical engineer and I have been using the FEM regularly for several years. However, I feel that I have never really understand the basis of the principle. I am using the heat equation to try to understand how the method is developed.
$$ \Delta{T} = q \\ T(\partial \Omega) = g $$
Firsts of all, I need to convert this PDE to the weak form, for which I multiply both sides by a test function $v$ and integrate. So first question, what are the conditions that $v$ must have? In some places I read "any test function", but in other places I read "compatible test functions".
$$ \int_\Omega {v\Delta Td\Omega} = \int_\Omega {vg d\Omega} \\ \int_\Omega v\Delta Td\Omega = -\int_\Omega {\nabla{v} \nabla{T}} d\Omega + \int_{\partial \Omega} v (\nabla T \cdot n) \partial \Omega = \int_\Omega vg d\Omega $$
We can chose $T = \sum_i^N T_i \cdot \phi_i$ and I chose j test functions $v_j = \psi_j$ to get j equations
$$ \sum_i^N T_i \int_\Omega {\nabla{\psi_j} \nabla{\phi_i}} d\Omega = \sum_i^N T_i \int_{\partial \Omega} \psi_j (\nabla \phi_i \cdot n) \partial \Omega - \int_\Omega \psi_j g d\Omega $$
In nearly all the articles I have read, they now assume homogeneous boundary conditions, and they chose shape functions $\phi_i = \psi_j$ vanishing at the boundary. With that selection, the boundary integral dissappears, and the boundary condition is enforced by construction, as the field $T$ has a value of 0 by construction on the boundary.
However, what happens when we have non-homogeneous boundary conditions? I guess we can select different shape functions for the field interpolation and the test functions ($\phi_i \neq \psi_j$), with $\psi_j$ vanishing at the boundary. That way, the boundary integral dissapears, but I still get equations for the values $T_i$ at the boundary, which I can impose while solving the system of equations.
However, this approach allows me to have "nodal" boundary conditions, but I am not sure how to proceed with "distributed" boundary conditions. How I could treat a region in my domain with a specified T value? Do I need to do the integral on the boundary?
Partial answer without going into too much details, assuming that you are content with the choice of test functions $v$ that are zero on the boundary. We will keep the choice $\phi_i = \psi_j$. The trick to deal with boundary condition is to split $T$ to two parts $$ T = T_0 + T_g $$ such that $T_0(\partial\Omega) = 0$ and $T_g(\partial\Omega) = g$ and find the solution for. We plug this into the so called weak formulation $$ \int_\Omega {\nabla{v} \nabla{T}} d\Omega = -\int_\Omega vg d\Omega + \int_{\partial \Omega} v (\nabla T \cdot n) \partial \Omega $$ and get $$ \int_\Omega {\nabla{v} \nabla{T_0}} d\Omega = -\int_\Omega vg d\Omega + \underbrace{\int_{\partial \Omega} v (\nabla T \cdot n) \partial \Omega}_{=0} - \int_\Omega {\nabla{v} \nabla{T_g}} d\Omega . $$ The second term is zero thanks to the choice of $v$ vanishing at the boundary. We are now almost at the same position as in the case of homogeneous boundary, only this time not for $T$ but for $T_0$.
The extra term with $T_g$ is a bit problematic in the sense that we do not know exactly $T_g$ looks like. The key observation is that it actually does not matter how it looks like as long as it satisfies $T_g(\partial\Omega) = g$. The convenient choice for implementation is $$ T_g := \sum_i \tilde g_i \phi_i^{boundary} $$ where $\phi_i^{boundary}$ are the basis functions that have nonzero values on boundary and $\tilde g_i$ are computed from $g(x)$ in a way that satisfies $T_g(\partial\Omega) = g$, up to some discretization error.
This might be a bit clearer if you write it down in terms of vectors and matrices you get after discretization. You effectively get a problem $$ A x = b $$ where $A$ is the stiffness matrix, $b$ is the volume source term and $x$ the unknown vector. We know that $x$ should have some part with prescribed values that come from the boundary condition, so we split $x$ in to two parts $$ x = x_0 + x_g $$ where $x_g$ comes from degrees of freedom on boundary and we know its values from the boundary condition and $x_0$ is the rest. We just plug in the split, get $$ Ax_0 = b - Ax_g $$ and solve for $x_0$.
If a region in your domain would have some prescribed values, you could use the same approach of splitting $T = T_0 + T_p$ where $T_p$ is the part with prescribed values and $T_0$ is the rest.