{find->get}_dofs in examples

[gdmcbain]

Following merger of #775 as discussed in #764.

[gdmcbain]

Comments on style welcome. Setting out a nice idiom is an important part of docs/examples; I'm trying to do the right thing here rather than push my own way of doing things.

[gdmcbain]

What's the right way to get all the boundary dofs as in

scikit-fem/docs/examples/ex13.py

Lines 43 to 47 in dea7871

	boundary_dofs = basis.find_dofs()
	u = basis.zeros()
	u[boundary_dofs['positive'].all()] = 1.
	u = solve(*condense(A, x=u, D=boundary_dofs))

Simple replacement of find with get looks wrong.

I do get the correct answer with

boundary_dofs = np.concatenate([basis.get_dofs(b) for b in mesh.boundaries])

but is that required?

[kinnala]

What's the right way to get all the boundary dofs as in

scikit-fem/docs/examples/ex13.py

Lines 43 to 47 in dea7871

	boundary_dofs = basis.find_dofs()
	u = basis.zeros()
	u[boundary_dofs['positive'].all()] = 1.
	u = solve(*condense(A, x=u, D=boundary_dofs))

Simple replacement of find with get looks wrong.

I do get the correct answer with

boundary_dofs = np.concatenate([basis.get_dofs(b) for b in mesh.boundaries])

but is that required?

I interpret "all DOFs" to mean the set of all DOFs on the boundary of the domain which is basis.get_dofs(). However, what you want to do is to get "all DOFs on all named boundaries". Could something explicit like

boundary_dofs = basis.get_dofs("ground") | basis.get_dofs("positive")

work? Another option would be to add support for

boundary_dofs = basis.get_dofs(["ground", "positive"])

as suggested in #764.

[kinnala]

Added support for get_dofs(('left', 'right')) in #778. Maybe that could be useful in some of these?

[gdmcbain]

Right, i do mean union rather than concatenation. Either proposed syntax is good (union of get-dofs of str or get-dof of collection of str). I'll try that.

[gdmcbain]

O. K., thanks, yes, #778 is nice. i'll use that. Ta.

[kinnala]

I'll fix the bugs in #778 asap.

[kinnala]

Now it should be fixed but I'll wait for the tests to pass.

[gdmcbain]

This just leaves ex36.

scikit-fem/docs/examples/ex36.py

Lines 163 to 180 in 9ee89ba

	ddofs = [
	basis["u"].find_dofs(
	{"left": mesh.facets_satisfying(lambda x: x[0] < 1.e-6)},
	skip=["u^2", "u^3"]
	),
	basis["u"].find_dofs(
	{"bottom": mesh.facets_satisfying(lambda x: x[1] < 1.e-6)},
	skip=["u^1", "u^3"]
	),
	basis["u"].find_dofs(
	{"back": mesh.facets_satisfying(lambda x: x[2] < 1.e-6)},
	skip=["u^1", "u^2"]
	),
	basis["u"].find_dofs(
	{"front": mesh.facets_satisfying(lambda x: np.abs(x[2] - 1.) < 1e-6)},
	skip=["u^1", "u^2"]
	)
	]

[kinnala]

Alright, the last one is slightly complicated, let's see if I can do something about it.

[kinnala]

Here is something:

modified   docs/examples/ex36.py
@@ -144,7 +144,12 @@ def volume(w):
     return J
 
 
-mesh = MeshTet().refined(2)
+mesh = MeshTet().refined(2).with_boundaries({
+    'left': lambda x: x[0] < 1.e-6,
+    'bottom': lambda x: x[1] < 1.e-6,
+    'back': lambda x: x[2] < 1.e-6,
+    'front': lambda x: np.abs(x[2] - 1.) < 1e-6,
+})
 uelem = ElementVectorH1(ElementTetP2())
 pelem = ElementTetP1()
 elems = {
@@ -160,36 +165,17 @@ du = basis["u"].zeros()
 dp = basis["p"].zeros()
 stretch_ = 1.
 
-ddofs = [
-    basis["u"].find_dofs(
-        {"left": mesh.facets_satisfying(lambda x: x[0] < 1.e-6)},
-        skip=["u^2", "u^3"]
-    ),
-    basis["u"].find_dofs(
-        {"bottom": mesh.facets_satisfying(lambda x: x[1] < 1.e-6)},
-        skip=["u^1", "u^3"]
-    ),
-    basis["u"].find_dofs(
-        {"back": mesh.facets_satisfying(lambda x: x[2] < 1.e-6)},
-        skip=["u^1", "u^2"]
-    ),
-    basis["u"].find_dofs(
-        {"front": mesh.facets_satisfying(lambda x: np.abs(x[2] - 1.) < 1e-6)},
-        skip=["u^1", "u^2"]
-    )
+D = [
+    basis["u"].get_dofs('front').all('u^3'),
+    basis["u"].get_dofs('left').all('u^1'),
+    basis["u"].get_dofs('bottom').all('u^2'),
+    basis["u"].get_dofs('back').all('u^3'),
 ]
 
-dofs = {}
-for dof in ddofs:
-    dofs.update(dof)
-
-du[dofs["left"].all()] = 0.
-du[dofs["bottom"].all()] = 0.
-du[dofs["back"].all()] = 0.
-du[dofs["front"].all()] = stretch_
+du[D[0]] = stretch_
 
 I = np.hstack((
-    basis["u"].complement_dofs(dofs),
+    basis["u"].complement_dofs(np.hstack(D)),
     basis["u"].N + np.arange(basis["p"].N)
 ))

[gdmcbain]

I'm pretty sure we can get rid of the x[0] < 1.e-6 stuff, that's just an idiom common in FEniCS but not at all a good idea for reasons argued above. (Meanwhile, I've launched #791 for predefining boundaries in the built-in mesh-generators.)