Problem in defining my exact solution

variational formulation
1

Hi,

I am trying to test my code for this aim I am using the method of manufactured solution as follows:

# the exact solution

Ue_x = Expression('Coef_u0*(1 + pow(sin(x[0]/R), 2) * pow(sin(x[1]/R), 2))', Coef_u0 = Coef_u0,
                  R = R, degree = basedegree + degreeCG)    # x-component of displacment
Ue_y = Expression('Coef_u0*(1 + pow(sin(x[0]/R), 2) * pow(sin(x[1]/R), 2))', Coef_u0 = Coef_u0,
                  R = R, degree = basedegree + degreeCG)    # y-component of displacment

Ke_1 = Expression(("(2*Coef_u0/R)*sin(x[0]/R)*sin(x[0]/R)*pow(sin(x[1]/R),2)","(2*Coef_u0/R)*sin(x[1]/R)*sin(x[1]/R)*pow(sin(x[0]/R),2)"), R = R, Coef_u0 = Coef_u0,
                  degree = basedegree + degreeCurl)  # the first row of displacement gradient
Ke_2 = Expression(("(2*Coef_u0/R)*sin(x[0]/R)*sin(x[0]/R)*pow(sin(x[1]/R),2)","(2*Coef_u0/R)*sin(x[1]/R)*sin(x[1]/R)*pow(sin(x[0]/R),2)"), R = R, Coef_u0 = Coef_u0,
                  degree = basedegree + degreeCurl)  # the second row of displacement gradient

def Stress_ES(Ke_1,Ke_2):
    def_grad = as_tensor([[1.0 + Ke_1[0], Ke_1[1]], [Ke_2[0], 1.0 + Ke_2[1]]])     # deformation gradient
    return mu*def_grad + (2*lam*ln(det(def_grad.T*def_grad)) - mu)*inv(def_grad.T)

Pe_1 = Stress_ES(Ke_1,Ke_2)[0,:]  # 1st row of stress
Pe_2 = Stress_ES(Ke_1,Ke_2)[1,:]  # 2nd row of stress
PP = Stress_ES(Ke_1,Ke_2)

def Body_f(PP):
    return -div(PP)

Be_x = Body_f(PP)[0]   # x-component of the body force for the exact solution
Be_y = Body_f(PP)[1]   # y-component of the body force for the exact solution

This is a 2D problem which I am trying to calculate the error regarding Displacement (U), Displacement Gradient (K) and Stress (P ) but I am getting this error:

ufl.log.UFLException: Cannot determine geometric dimension from expression.

The source of error which comes to my mind might be:

  1. The way I have defined Exact Solution for my stress:

Pe_1 and Pe_2

  1. The way I am extracting row one and row two of my body force:
Be_x = Body_f(PP)[0] and  Be_y = Body_f(PP)[1]

Thank you in advance.

2

Consider using SpatialCoordinates instead.
I’ve added an example with dummy parameters below:

from dolfin import *
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 1)

x = SpatialCoordinate(mesh)
Coef_u0 = Constant(5)
R = 3
lam = 2
Ue_x = Coef_u0*(1 + pow(sin(x[0]/R), 2) * pow(sin(x[1]/R), 2))   # x-component of displacment
Ue_y = Coef_u0*(1 + pow(sin(x[0]/R), 2) * pow(sin(x[1]/R), 2))
Ke_1 = as_vector(((2*Coef_u0/R)*sin(x[0]/R)*sin(x[0]/R)*pow(sin(x[1]/R),2),
                  (2*Coef_u0/R)*sin(x[1]/R)*sin(x[1]/R)*pow(sin(x[0]/R),2)))
Ke_2 = as_vector(((2*Coef_u0/R)*sin(x[0]/R)*sin(x[0]/R)*pow(sin(x[1]/R),2),
                  (2*Coef_u0/R)*sin(x[1]/R)*sin(x[1]/R)*pow(sin(x[0]/R),2)))
def_grad = as_tensor([[1.0 + Ke_1[0], Ke_1[1]], [Ke_2[0], 1.0 + Ke_2[1]]])     # deformation gradient
mu = 0.01
Stress_ES = mu*def_grad + (2*lam*ln(det(def_grad.T*def_grad)) - mu)*inv(def_grad.T)

Pe_1 = Stress_ES[0,:]  # 1st row of stress
Pe_2 = Stress_ES[1,:]  # 2nd row of stress
PP = Stress_ES

def Body_f(PP):
    return -div(PP)

Be_x = Body_f(PP)[0]   # x-component of the body force for the exact solution
Be_y = Body_f(PP)[1]   # y-component of the body force for the exact solution
print(Be_x)
1 Like
3

Hi Jorgen,

Thanks for your response.
I’ve modified my script and saw this error there.

Reason: Cannot be created from subspace. Consider collapsing the function space.
*** Where: This error was encountered inside Function.cpp.

This is the way I am imposing my BCs:

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

#      Mark bottom boundary facets as subdomain 0
class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-12   # tolerance for coordinate comparisons
        return on_boundary and abs(x[1]) < tol
Gamma_B = BottomBoundary()
Gamma_B.mark(boundary_parts, 0)

#       Mark top boundary facets as subdomain 1
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-12   # tolerance for coordinate comparisons
        return on_boundary and abs(x[1] - 1) < tol
Gamma_T = TopBoundary()
Gamma_T.mark(boundary_parts, 1)

#       Mark left boundary as subdomain 2
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-12   # tolerance for coordinate comparisons
        return on_boundary and abs(x[0]) < tol
Gamma_L = LeftBoundary()
Gamma_L.mark(boundary_parts, 2)

#       Mark right boundary as subdomain 3
class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-12   # tolerance for coordinate comparisons
        return on_boundary and abs(x[0] - 1) < tol
Gamma_R = RightBoundary()
Gamma_R.mark(boundary_parts, 3)

#       Bottom boundary function
def bottom_f(x, on_boundary):
    tol = 1E-12   # tolerance for coordinate comparisons
    return on_boundary and abs(x[1]) < tol


#       BCs for the initial guess and the increment in each Newton's iteration

bcs = [DirichletBC(Z.sub(0), Ue_x, bottom_f),
       DirichletBC(Z.sub(1), Ue_y, bottom_f)]

Then I modified to : (Using “collapse”)

#BCs for the initial guess and the increment in each Newton’s iteration
bcs = [DirichletBC(Z.sub(0).collapse(), Ue_x, bottom_f),
DirichletBC(Z.sub(1).collapse(), Ue_y, bottom_f)]

and see this error again:

Error: Unable to define linear variational problem a(u, v) = L(v) for all v.
*** Reason: Expecting the boundary conditions to live on (a subspace of) the trial space.

4

I suggest to interpolate Ue_x and Ue_y into the appropriate collapsed function spaces.