node ordering and include assumptions to avoid occurences of sympy Piecewise

[stefanhiemer]

Hi there,

two questions:

i) Are there any utility functions to change the default node ordering to some of the common FE-codes? If no and you have ideas what the proper ordering algorithms are, I would be happy to contribute.

ii) I try to derive the matrix for a monomial of order n. For some predetermined n, no problem, but I wanted to do it in a general way. The integral looks like this N.T @ (N.T @ c)^(n-1) N.T where N are the shape functions and c is the vector of nodal variables of interest (e. g. concentration). As a result I get a very ugly Piecewise defined function:

Matrix([
[Piecewise((2*l1*(n**2*u1**n*u2**2 - 2*n**2*u1**(n + 1)*u2 + n**2*u1**(n + 2) + 3*n*u1**n*u2**2 - 4*n*u1**(n + 1)*u2 + n*u1**(n + 2) + 2*u1**n*u2**2 - 2*u2**(n + 2))/(n*(n**2*u1**3 - 3*n**2*u1**2*u2 + 3*n**2*u1*u2**2 - n**2*u2**3 + 3*n*u1**3 - 9*n*u1**2*u2 + 9*n*u1*u2**2 - 3*n*u2**3 + 2*u1**3 - 6*u1**2*u2 + 6*u1*u2**2 - 2*u2**3)), Ne(u1, u2)), (2*l1*u2**(n - 1)/3, True)),                                             Piecewise((2*l1*(n*u1*u2**(n + 1) - n*u1**(n + 1)*u2 + n*u1**(n + 2) - n*u2**(n + 2) + 2*u1*u2**(n + 1) - 2*u1**(n + 1)*u2)/(n*(n**2*u1**3 - 3*n**2*u1**2*u2 + 3*n**2*u1*u2**2 - n**2*u2**3 + 3*n*u1**3 - 9*n*u1**2*u2 + 9*n*u1*u2**2 - 3*n*u2**3 + 2*u1**3 - 6*u1**2*u2 + 6*u1*u2**2 - 2*u2**3)), Ne(u1, u2)), (l1*u2**(n - 1)/3, True))],
[                                           Piecewise((2*l1*(n*u1*u2**(n + 1) - n*u1**(n + 1)*u2 + n*u1**(n + 2) - n*u2**(n + 2) + 2*u1*u2**(n + 1) - 2*u1**(n + 1)*u2)/(n*(n**2*u1**3 - 3*n**2*u1**2*u2 + 3*n**2*u1*u2**2 - n**2*u2**3 + 3*n*u1**3 - 9*n*u1**2*u2 + 9*n*u1*u2**2 - 3*n*u2**3 + 2*u1**3 - 6*u1**2*u2 + 6*u1*u2**2 - 2*u2**3)), Ne(u1, u2)), (l1*u2**(n - 1)/3, True)), Piecewise((2*l1*(-n**2*u1**2*u2**n + 2*n**2*u1*u2**(n + 1) - n**2*u2**(n + 2) - 3*n*u1**2*u2**n + 4*n*u1*u2**(n + 1) - n*u2**(n + 2) - 2*u1**2*u2**n + 2*u1**(n + 2))/(n*(n**2*u1**3 - 3*n**2*u1**2*u2 + 3*n**2*u1*u2**2 - n**2*u2**3 + 3*n*u1**3 - 9*n*u1**2*u2 + 9*n*u1*u2**2 - 3*n*u2**3 + 2*u1**3 - 6*u1**2*u2 + 6*u1*u2**2 - 2*u2**3)), Ne(u1, u2)), (2*l1*u2**(n - 1)/3, True))]])

(I guess) sympy tries to deal with the case n =< 0, but this case is of no interest (at least to me). My question is now how can I switch of this behaviour as I am only interested in the case of n being an integer and greater/equal 2. Below is the full code to check.

mononomial.txt

Thanks and Regards
Stefan

[mscroggs]

(i) We don't have any functions to permute the order of DOFs to match libraries (although the ordering we use matches FEniCSx/Basix), but I'd be happy to have these added. Functions called something symfem.permuting.LIBRARY_NAME that take an element as input then return a permutation might be a way to do this. Which library/libraries do you have in mind?

(ii) You can use assumptions in Sympy, although it's a while since I used them. Looks like something like sympy.Symbol("n", positive=True, integer=True) might be what you want. There's more examples at https://docs.sympy.org/latest/guides/assumptions.html.

[stefanhiemer]

Sorry for my late answer as this actually took me down a rabbit hole.

i) I currently work with ElmerFEM

https://en.wikipedia.org/wiki/Elmer_FEM_solver
https://www.elmerfem.org/blog/

and would try to be compatible with Abaqus (as much as I can without a license) although I openly admit I have close to no experience with the latter. So as soon as I have something I would get back at you regarding this point.

ii) If you check the txt-file that is already what I did (line 40), but I found out what the problem is: The assumption n>=2 does not trigger the Piecewise definition, but the nodal variables (symbols u1 and u2) are the culprits as the integral has several different expression/simplifications when assuming them to be related to the element dimensions. This cannot be suppressed not matter what assumptions you apply as apparently during the expression evaluation/integration sympy does not (yet) care about the assumption context. Sympy is in the process of switching between the old assumption system linked to individual symbols which is taken care of during evaluation e. g. sympy.Symbol("n", positive=True, integer=True). The new system with the assumption context e. g.

with assuming(Q.ge(n, 2), ): integrand = integrand.integral(ref,x)
is not yet considered during expression evaluation as far as I can see. So no matter which inequalities I apply in the assuming expression, the resulting integrand expression is always the same. So as far as I see, nothing can be done at this point in time as this appears to be a sympy problem.

[stefanhiemer]

Hi there

I’ve started working on an automatic abaqus-like node ordering and have two routines:

i) a counter-clockwise sorting that orders nodes by angle/radius (and z in 3D). For the first routine to be finished I must first be classified by topological entity (vertex / edge / face / cell interior) before any geometric ordering. Does Symfem already expose information about which basis functions / DOFs are associated with which entity (vertex, edge, face, cell)?

https://github.com/stefanhiemer/topoptlab/blob/feature/inverse-homogenization/topoptlab/symbolic/utils.py

ii) for finding the nodal position of the basis functions I use (solve(phi_i - 1 == 0)) to extract nodal coordinates by using the delta property:

https://github.com/stefanhiemer/topoptlab/blob/feature/inverse-homogenization/topoptlab/symbolic/cell.py

However this is rather painful as one has to add a number of inequalities, exlude roots, etc. pp. Is there any existing Symfem functionality that directly returns nodal coordinates, or is solving the delta property the current approach to take?

Thanks and regards
Stefan