What's meaning of the error message, "The provided data did not satisfy the prerequisites"?

5 months ago by
I have two questions: 1.why I get a NaN? 
2. After uncommenting, What's meaning of the error message, "The provided data did not satisfy the prerequisites"?

When comment 
result1 = mysolver_1(mesh, f1)​

, I get the result and get a NaN. Or, I get the below error:
*** Error: Unable to compute matrix factorisation.
*** Reason: The provided data did not satisfy the prerequisites.
*** Where: This error was encountered inside EigenLUSolver.cpp.
*** Process: 0
*** DOLFIN version: 2017.1.0
*** Git changeset: unknown

Here is my simple source code. Forgive me it is also long, but you don't need to understand the whole code.

import pdb
from fenics import *
from mshr import *
import numpy as np
from pprint import pprint
from scipy import sparse

# form source terms
def all_Expression():
    cell = 'tetrahedron'
    import sympy as sp
    x, y, z = sp.symbols('x[0] x[1] x[2]')
    lamda = 1.0
    p = x - 0.5
    U = (x*(x-1) * y*(y-1) * z*(z-1) * (y-0.5)*(z-0.5))**3
    p_x = p.diff(x, 1)
    p_y = p.diff(y, 1)
    p_z = p.diff(z, 1)
    U_x = U.diff(x, 1)
    U_y = U.diff(y, 1)
    U_z = U.diff(z, 1)
    Laplace_U = sum(U.diff(xi, 2) for xi in (x, y, z))
    Laplace_Laplace_U = sum(Laplace_U.diff(xi, 2) for xi in (x, y, z))
    #---------------source terms-----------------------
    p_e = Expression(sp.printing.ccode(p), degree=1, cell=cell)
    f1 = Expression((sp.printing.ccode(- p_x), \
                    sp.printing.ccode(- p_y), \
                    sp.printing.ccode(- p_z)), degree=1, cell=cell)
    f2 = Expression(sp.printing.ccode(- Laplace_Laplace_U - lamda - p), degree=4, cell=cell)

    return p_e, f1, f2

def mysolver_1(mesh, f):
    parameters.linear_algebra_backend = 'Eigen'
    h = mesh.hmax()

    DG = FunctionSpace(mesh, 'DG', 0)
    RT = FunctionSpace(mesh, 'RT', 1)
    p = TrialFunction(DG)
    v = TestFunction(RT)
    bc = DirichletBC(RT, Constant([0.0, 0.0, 0.0]), DomainBoundary())
    a = p * div(v) * dx
    L = inner(f, v) *dx
    a_ = p * dx

    A, b = assemble_system(a, L, bc)
    A_ = assemble(a_)
    rows, columns, values = as_backend_type(A).data()
    csrA = sparse.csr_matrix((values, columns, rows))
    AA = sparse.vstack([csrA, A_.array().T])
    b_ = list(b.array()); b_.append(0)
    bb = np.array(b_)

    from scipy.sparse.linalg import lsqr
    pp = lsqr(AA, bb)[0]
    p = Function(DG)

    plot(p, title='solution p')
    return p, h

def mysolver_2(mesh, f, p):
    h = mesh.hmax()

    CG1 = FiniteElement('CG', mesh.ufl_cell(), 1)
    R = FiniteElement('R', mesh.ufl_cell(), 0)
    element = MixedElement([CG1, R])
    V = FunctionSpace(mesh, element)
    bc = DirichletBC(V.sub(0), Constant(0.0), DomainBoundary())

    phi, lamda = TrialFunctions(V)
    psi, mu = TestFunctions(V)

    a = inner(grad(phi), grad(psi))*dx - lamda*psi*dx \
            - phi*mu*dx
    L = p*psi*dx + f*psi*dx

    U = Function(V)
    solve(a == L, U, bc)

    phi, lamda = U.split()
    plot(phi, title='phi')
    plot(lamda, title='lamda')
    return phi, lamda, h

if __name__ == '__main__':
    p_e, f1, f2 = all_Expression()
    domain = Box(Point(0,0,0), Point(1,1,1)) - Box(Point(0,0,0), Point(1,0.5,0.5))
    mesh = generate_mesh(domain, 1)

    for i in range(3):
        #----------------------------------------error here----------------------------------------------------
        result1 = mysolver_1(mesh, f1)

        result2 = mysolver_2(mesh, f2, p_e)
        mesh = refine(mesh)


Community: FEniCS Project
See https://www.allanswered.com/post/emre/how-to-ask-a-good-question-and-get-an-answer/ for advice on asking questions.
written 5 months ago by Garth Wells  
I am sorry for the inappropriate question.
Now I make the question clear.
Thanks a lot.
written 5 months ago by fanronghong  

1 Answer

5 months ago by

The problem actually occurs for me in solver 2, rather than solver 1.

I believe that the reason is that the mesh you generate with "generate_mesh(domain, 1)" is so coarse that every node is on the boundary, so they are all constrained by DirichletBC. This results in a locking for the CG1-R mixed space and thus a singular system.


Please login to add an answer/comment or follow this question.

Similar posts:
Search »