Dear all, In the code attached, the objective is to copy the surface values of the 3D function phi_s, to the interior nodes of the 3D function phi_3s such that phi_3s is equal to the top values of phi_2s at each z. For this, I want to solve dzz(phi) = 0, with boundary conditions : phi_3s = phi_s(top) at the bottom z=0 and top z=H0. solving the equation: phi.dx(2)*v.dx(2) = 0 (where phi is the trial function and v the test function) did not work, so instead I solve: grad(phi).grad(v) = phi_3s.dx(0)v.dx(0) + phi_3s.dx(1) v.dx(1), with phi_s as boundary condition. I have two questions concerning this code: 1) First, when I change CG1 to CG2, then the value of phi_3s is not exactly the same as the one of phi_s at the surface. Where does this difference come from? 2) Is there a way to define a boundary condition at the bottom of the domain V (at z=0), using the value of a function at the top of this domain? Basically, I would like to use the top surface value of phi_s as a boundary condition both at the bottom z=0 and at the top z=h0 of the domain. Is that possible? Thanks and best regards, Floriane