Matrix valued coefficient function

More
3 years 10 months ago - 3 years 10 months ago #3474 by Ryanleaf
Hello everybody,

Sorry, I am fairly new to python and NGSolve, I think I have a relatively simple question.
I am trying to compute the RHS of the surface stokes equations by a given analytic solution.
The geometry I have is just a simple unit sphere. So I wrote the following operators to compute RHS of surface stokes equations:
Code:
def Coef_Grad(v): return CoefficientFunction(tuple([v[i].Diff(w) for w in [x, y, z] for i in [0, 1, 2]]), dims=(d, d)).trans def Coef_Matrix_Multi_Matrix(P,Q): Q_list=[[Q[0,0], Q[0,1], Q[0,2]], [Q[1,0], Q[1,1], Q[1,2]], [Q[2,0], Q[2,1], Q[2,2]]] P_list=[[P[0,0], P[0,1], P[0,2]], [P[1,0], P[1,1], P[1,2]], [P[2,0], P[2,1], P[2,2]]] return CoefficientFunction(tuple([sum([P_list[i][k]*Q_list[k][j] for k in range(d)]) for j in [0, 1, 2] for i in [0, 1, 2]]), dims=(d, d)) def Coef_Grad_Gamma(v): # set projection matrix Proj = [[0.5*(1-x**2), -x*y, -x*z], [0, 0.5*(1-y**2), -y*z], [0, 0, 0.5*(1-z**2)]] proj = CoefficientFunction(tuple([Proj[w][i] for w in [0, 1, 2] for i in [0, 1, 2]]), dims=(d, d)) proj = proj+proj.trans v_coef = CoefficientFunction(tuple([v[i] for i in [0, 1, 2]])) return Coef_Matrix_Multi_Matrix(Coef_Matrix_Multi_Matrix(proj, Coef_Grad(v_coef)), proj) def Coef_Eps_Gamma(v): return 0.5 * (Coef_Grad_Gamma(v) + Coef_Grad_Gamma(v).trans) def Coef_Div_Gamma(v): v_list = [[v[0,0], v[0,1], v[0,2]], [v[1,0], v[1,1], v[1,2]], [v[2,0], v[2,1], v[2,2]]] return CoefficientFunction(tuple([Trace(Coef_Grad_Gamma(v_list[i])) for i in range(d)]))
I didn't find how to modified the matrix valued coefficient function after define it, sorry I used the list.
The problem I have now is When I tried to compute
Code:
coef_b = Coef_Div_Gamma(-2 * mu[2] * Coef_Eps_Gamma(vel_exact[2]))
And assemble corresponding RHS in the algebraic system
Code:
f += SymbolicLFI(lset_if, form=Pmat*coef_b * v[2])
My code became very very slow on the assembling step. I think the main reason might be (Coef_Div_Gamma(v)) operator, but I don't know what is wrong with this part. I think it took days to assemble, any reasons caused this?
Where, Pmat defined as,
Code:
Pmat = Id(3) - OuterProduct(n_lset,n_lset)
v[2] belongs to
Code:
Vhbase = VectorH1(mesh, order=order, dirichlet="bottom|top|left|right|front|back") VhG=Compress(Vhbase,GetDofsOfElements(Vhbase, ci.GetElementsOfType(IF)))
thanks a lot,
Ryan
Last edit: 3 years 10 months ago by Ryanleaf.
More
3 years 10 months ago - 3 years 10 months ago #3475 by Ryanleaf
Sorry problem solved by reading one of the previous posts.
Problem solved by directly using matrix-valued coefficient function multiplication.
Last edit: 3 years 10 months ago by Ryanleaf.
Time to create page: 0.098 seconds