# Have a matrix solution from affine operator

**URL:** <https://discourse.julialang.org/t/have-a-matrix-solution-from-affine-operator/100026>\
**Category:** New to Julia\
**Tags:** question, gridap\
**Created:** [June 8, 2023, 12:14am UTC](https://discourse.julialang.org/t/have-a-matrix-solution-from-affine-operator/100026 "2023-06-08T00:14:11Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![mary](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mary/32/53139_2.png) [@mary](https://discourse.julialang.org/u/mary)\
**Post date:** [June 8, 2023, 12:14am UTC](https://discourse.julialang.org/t/have-a-matrix-solution-from-affine-operator/100026/1 "2023-06-08T00:14:11Z")

</div>

Hi,  
I am solving a linear elasticity problem using affine and solve and the output of this function is single field function. I need the matrix form of the solution to do some calculation on it. can you tell me how can i do this? i write it tofind the answer in matrix format but it gives me zero which is not correct. even if you can modify the matrix code instead to get right answers i would appreciate your help.  
the lines using affine are as below

```julia
function stepDisp(fem_params,pth)
    uApp1(x) = VectorValue(0.0,0.0)
    uApp2(x) = VectorValue(0.0,0.0)
    uApp3(x) = VectorValue(0.0,-vAppMax)
    U_Disp = TrialFESpace(V0_Disp ,[uApp1 ,uApp2 ,uApp3])
    A_Disp(u,v,pth) = ((p->Em(p))∘pth) * (ε(v) ⊙ (σfun∘(ε(u))))
    a_Disp(u,v) =∫(A_Disp(u,v,pth))fem_params.dΩ
    b_Disp(v) = 0.0
    op_Disp = AffineFEOperator(a_Disp ,b_Disp ,fem_params.U_Disp ,fem_params.V0_Disp)
    uh_out = solve(op_Disp)
    return uh_out
end
vApp = vAppMax
p0 = ones(fem_params.np)
pf_vec = Filter(p0;r, fem_params)
pfh = FEFunction(fem_params.Pf, pf_vec)
pth = (pf -> Threshold(pf; β, η)) ∘ pfh
A_mat = MatrixA(pth; fem_params)
u_vec = stepDisp(fem_params,pth)

```

and the matrix format is like below

```julia
vApp = vAppMax
p0 = ones(fem_params.np)
pf_vec = Filter(p0;r, fem_params)
pfh = FEFunction(fem_params.Pf, pf_vec)
pth = (pf -> Threshold(pf; β, η)) ∘ pfh
A_mat = MatrixA(pth; fem_params)
b_vec = assemble_vector(0.0, fem_params.V0_Disp)
u_vec = A_mat \ b_vec

```

and i want to calculate the lines below after solving finite element that i have provided.

```julia
function MatrixOf(fem_params)

    return assemble_matrix(fem_params.U_Disp, fem_params.V0_Disp) do u, v
              0.5*(∫((∇(u))' ⊙ (C_mat ⊙ ∇(v)))fem_params.dΩ)
    end
end
O_mat = MatrixOf(fem_params)
    u_vec' * O_mat * u_vec

```

and this gives me assertation error of number of arrays and triangulation.  
thank you all