Find true degree of freedom

[BBBBBbruce]

Hello, I am new to computational mechanics, so hope my question is not too silly.

I am trying to build a simple direct solver, instead of using the iterative solver. The reason is that the direct solver is usually faster than the iterative solver.

I am using ex2 as an example. My code works ok when I disable the mesh refinement(section 4), I compared it with the result calculated by PCG, they are quite similar.

But when I enable the mesh refinement, my _A(defined below, the reduced stiffness matrix) becomes singular. I look the code of GetEssentialTrueDofs(), but I could grasp what exactly is done to calculate the true dofs? Also, what boundary conditions should be added to prevent the singular matrix problem?

Here is my code of direct solver:

    //================ Direct solver ==========================
    //input A, B, ess_tdof_list;
    //output __X

    //reducing B to _B
    Vector _B;
    _B.SetSize(B.Size()-ess_tdof_list.Size());
    int j = 0;
    ess_tdof_list.Print();
    for (int i = 0; i < _B.Size(); i++) {
        if (i+j == ess_tdof_list[j]) {
            j++;
        }
        _B[i] = B[i + j];
    }
    
    //reducing A to _A
    DenseMatrix A_d = *A.ToDenseMatrix();
    DenseMatrix _A = DenseMatrix(A.Size() - ess_tdof_list.Size());
    
    int m = 0;
    int n = 0;

    for (int i = 0; i < _A.Size(); i++) {
        for (int k = 0; k < _A.Size(); k++) {
            if (i + m == ess_tdof_list[m] && k + n == ess_tdof_list[n]) {
                m++;
                n++;
            }
            else if (i + m == ess_tdof_list[m] && k + n != ess_tdof_list[n]) {
                m++;
            }
            else if (i + m != ess_tdof_list[m] && k + n == ess_tdof_list[n]) {
                n++;
            }

            _A.Elem(i, k) = A_d.Elem(i + m, k + n);
        }
        n = 0;
    }

    // initialise zero vector
    Vector _X;
    _X.SetSize(_B.Size());
    for (int i = 0; i < _B.Size(); i++)
        _X[i] = 0.0;
    
    _A.Invert();
    _A.AddMult( _B,_X);
    
    //Pad with fixed points;
    Vector __X;
    __X.SetSize(_X.Size() + ess_tdof_list.Size());
    j = 0;
    for (int i = 0; i < __X.Size(); i++) {
        if (i != ess_tdof_list[j]) {
            __X[i] = _X[i - j];

        }
        else {
            j++;
            __X[i] = 0.0;
        }
    }

[v-dobrev]

I see that you are using a direct dense matrix solver -- this will become very slow very quickly with refinement -- the inversion has complexity O(n^3) where n is the size of linear system.

If you want to use a direct solver, a much better solution is to use a direct sparse solver such as UMFPACK from the SuiteSparse library. We have support for UMFPACK in MFEM and, in fact, ex2.cpp will use that instead of PCG when it is build with SuiteSparse, see here

mfem/examples/ex2.cpp

Lines 202 to 206 in 1d03963

	// 11. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
	UMFPackSolver umf_solver;
	umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
	umf_solver.SetOperator(A);
	umf_solver.Mult(B, X);

[BBBBBbruce]

Hello, Thanks for the reply, I have get it to work. But I still have a theoretical problem: how to apply the true dofs?

From the teaching materials I read, the key is to "eliminate them". The naive way is to complete remove the corresponding rows and columns(called it reduced version) from the stiffness matrix. However, this linear system is not solvable in this case.

Additionally, from my investigation, mfem is set the true dofs to zeros(call it set-to-zero version). And the linear system is solvable.

This is where I got confused. In summary, my two questions are: 1. why mfem setting the true dofs to zero?(could you recommand any resources?) 2. I tested the condition number of the stiffness matrix in both reduced version and set-to-zero version. They are similar-- this means the invertibility are similar. Then, why is set-to-zero version is solvable?

Thank you for the help

[tzanio]

The way we impose BCs is briefly discusses here:

#1720 (comment)

Is this helpful?

[BBBBBbruce]

Thank you, it is a very detailed explanation. It is very helpful, but apologised I am very beginning of computational mechanics, I need a bit time to understand and see if it can solve my questions. Can I leave this open for a few days?

[BBBBBbruce]

The way we impose BCs is briefly discusses here:

#1720 (comment)

Is this helpful?

This answers my questions, Thanks!