search

kinnala / scikit-fem

Simple finite element assemblers
https://scikit-fem.readthedocs.io
BSD 3-Clause "New" or "Revised" License
500 stars 79 forks source link

get_dofs(nodes: Any) #1078

Closed JarosT55 closed 9 months ago

JarosT55 commented 9 months ago

Hi, I am trying to call the _getdofs() function on listed meshes boundary nodes like: _basis.get_dofs(nodes={'bnodes_a','bnodesb'}) The result is 'NotImplementedError()' The bnodes are not on a trivial path; thus, extracting them by position query would be difficult. I appreciate any suggestions. Regards MichaelT

kinnala commented 9 months ago

We don't have a way to tag (give a name) for a set of nodes, and nodes kwarg is for extracting DOFs for single point at a time.

Perhaps you mean: facets={'bnodes_a', 'bnodes_b'}) ?

gdmcbain commented 9 months ago

How is your basis defined?

Does ex28 work for you? It loads a mesh from a file, where the mesh has named boundaries.

https://github.com/kinnala/scikit-fem/blob/44916dc9551dcf5a3b12586080d61e18dc39dabe/docs/examples/ex28.py#L72

then accesses a subset of those

https://github.com/kinnala/scikit-fem/blob/44916dc9551dcf5a3b12586080d61e18dc39dabe/docs/examples/ex28.py#L107

gdmcbain commented 9 months ago

Oh I missed the bit about nodes. Do you actually want nodes? Or just to access predefined boundaries?

JarosT55 commented 9 months ago

Thank you All for the prompt response. I have not expected them ... so fast!!!! as per doc nodes kwarg is listed in the func description _getdofs(facets: Any | None = None, elements: Any | None = None, nodes: Any | None = None, skip: List[str] | None = None)

I will try playing with ex28.py

Regards MichaelT

kinnala commented 9 months ago

Perhaps you can try also: help(basis.get_dofs)

JarosT55 commented 9 months ago

Hi All,

@ _D = basis["heat"].getdofs([b for b in mesh.boundaries if b.endswith("-inlet")])

I use a similar syntax on the mesh. From a custom utility plot, I know it's every node, number, position, etc. I must admit that I am 2 weeks old skfem-infant still crawling on all four; thus, there are many mysteries at the front of my eyes. One confusion is that the list of nodes from: _mesh.boundaries['bnodesa'] produces a different list than the one from (for the ElementTriP1()) _basis.get_dofs(['bnodes_a']).nodalix ... but hey! _basis.get_dofs(['bnodes_a']).facetix ... is the perfect match!!! … and even better for TriMorley() I have got _‘u_n’_s _basis.get_dofs(['bnodes_a']).facet[‘un’]

Thus, ... there are a lot of puzzles to put into place still ... in the skfem environment, which looks accommodating and promising!

@ help(basis.get_dofs)

... is very useful!

Regards MichaelT

kinnala commented 9 months ago

In scikit-fem, you will need to carefully distinguish between different topological entities (e.g., nodes, facets) of the mesh AND "degrees-of-freedom" (DOF) which are, in general, different. Especially if you use an element such as ElementTriMorley which has both, nodal DOFs and facet DOFs.

For example, mesh.boundaries['bnodes_a'] does not produce a list of node indices but a list of facet indices. On the other hand, basis.get_dofs(['bnodes_a']).XXX produces a list of DOFs. Some discussion about this topic can be found from the documentation: https://scikit-fem.readthedocs.io/en/latest/advanced.html#indexing-of-the-degrees-of-freedom

JarosT55 commented 9 months ago

![Uploading ex28_basis_heat.png…]() Hi, As per gdmcbain advice, I have played a little bit with ex28. The extraction of boundaries has worked as a charm:

        _boundaryIdxes = basis.get_dofs([boundaryName]).nodal_ix     
        ax.plot(*mesh.p[:, boundaryIdxes], color=boundaryColor, linewidth=2*boundaryLine, )_

However, while attempting to visualize subdomains … I encountered some problems I could not comprehend. I assumed that each subdomain index addresses facet holding indexes of two nodes (defining the facet segment). I used the code like:

domainIdxes = mesh.subdomains[domainName] sIdxes = mesh.facets[0][domainIdxes] eIdxes = mesh.facets[1][domainIdxes] seIdxes = [i for indexes in zip(sIdxes, eIdxes) for i in indexes] ax.plot( *mesh.p[:, seIdxes], domainColor, markersize=domainMarker, label=domainName, linewidth=subdomainlinewidth )

The resulting plot … is a little bit messy! I do not know why. Is my assumption correct for organizing data? Is there a better/proper way of extracting subdomains? I am attaching ex28_basis_heat.png file with a visual representation of my test play.

Regards and HNY MichaelT ex28_basis_heat

kinnala commented 9 months ago

Subdomains are element indices, not facet indices. This does not make sense since you are slicing facets using element indices:

domainIdxes = mesh.subdomains[domainName] sIdxes = mesh.facets[0][domainIdxes]

JarosT55 commented 9 months ago

ooohhhhh ! Thanks

JarosT55 commented 9 months ago

uuuuuuFfffffff !!!

The code is ... well ... certainly not very efficient ... but it ![Uploading ex28_basis_heat.png…]() works.

_for subdomain, domainColor in zip(mesh.subdomains, colLst):

Elements centroids

domain_mesh = mesh.restrict(mesh.subdomains[subdomain]) domain_basis = Basis(domain_mesh, basis.elem) domain_basis.get_dofs(elements= lambda x: ax.plot(x[0], x[1], domainColor, markersize=markersizeM, label='elements'))

Elements facets

mesh_t2f = mesh.t2f[:, mesh.subdomains[domainName]]domain_facets = mesh.facets[:,mesh_t2f]
sIdxes = domain_facets[0] eIdxes = domainfacets[1] seIdxes = [i for indexes in zip(sIdxes, eIdxes) for i in indexes] ax.plot( *mesh.p[:, seIdxes], domainColor, markersize=domainMarker, label=domainName, linewidth=linewidthM )

Thank you for the valuable advices.

Regards MichaelT

ex28_basis_heat