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:

```
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

```
coef_b = Coef_Div_Gamma(-2 * mu[2] * Coef_Eps_Gamma(vel_exact[2]))
```

And assemble corresponding RHS in the algebraic system

`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,

`Pmat = Id(3) - OuterProduct(n_lset,n_lset)`

v[2] belongs to

```
Vhbase = VectorH1(mesh, order=order, dirichlet="bottom|top|left|right|front|back")
VhG=Compress(Vhbase,GetDofsOfElements(Vhbase, ci.GetElementsOfType(IF)))
```

thanks a lot,

Ryan