Assigning dirichlet BC to a component of a vector field

[bhaveshshrimali]

In a way this is related to #455 but may be relevant in general too. What's the best way to select dofs that belong to a specific component of the solution?

For instance, consider a unit cube where I want to specify u_1 = 0 on the left (x=0), u_2 = 0 on the bottom (y=0) and u_3 = 0 on the back (z=0) faces. Note that the remaining two components on each face are left unconstrained (physically, like in uniaxial tension). Currently I am tempted to think of it in a brute force way:

from numpy import hstack, unique
stretch_ = 1.5
dofs = {
    "left":basis["u"].get_dofs(lambda x: x[0] == 0.),
    "bottom":basis["u"].get_dofs(lambda x: x[1] == 0.),
    "back":basis["u"].get_dofs(lambda x: x[2] == 0.),
    "front":basis["u"].get_dofs(lambda x: x[2] == 1.)
}

du[dofs["left"].nodal["u^1"]] = 0.
du[dofs["left"].edge["u^1"]] = 0.

du[dofs["bottom"].nodal["u^2"]] = 0.
du[dofs["bottom"].edge["u^2"]] = 0.

du[dofs["back"].nodal["u^3"]] = 0.
du[dofs["back"].edge["u^3"]] = 0.

du[dofs["front"].nodal["u^3"]] = stretch_
du[dofs["front"].edge["u^3"]] = stretch_

I = basis["u"].complement_dofs(dofs) #this selects the complement of everything in `dofs`

I = hstack((
    I, dofs["left"].nodal["u^2"],
    dofs["left"].nodal["u^3"],
    dofs["left"].edge["u^2"],
    dofs["left"].edge["u^3"],
    dofs["bottom"].nodal["u^1"],
    dofs["bottom"].nodal["u^3"],
    dofs["bottom"].edge["u^1"],
    dofs["bottom"].edge["u^3"],
    dofs["back"].nodal["u^1"],
    dofs["back"].nodal["u^2"],
    dofs["back"].edge["u^1"],
    dofs["back"].edge["u^2"],
    dofs["front"].nodal["u^1"],
    dofs["front"].nodal["u^2"],
    dofs["front"].edge["u^1"],
    dofs["front"].edge["u^2"]
))

I = unique(I)

I am still thinking if the above is correct (as in it doesn't duplicate dofs or selects dofs that have dirichlet boundary specified).

Naturally, what would be the best way to assign dofs such that only the ones corresponding to a component are assigned and the complement_dofs method takes care of the rest ? (Is there any other method instead, that I should use?)

[gdmcbain]

In ex27 I had to do something similar: impose a Dirichlet condition on the velocity at the inlet. That used find_dofs rather than get_dofs as there was a period when the latter was slated for deprecation. Anyway, that is:

scikit-fem/docs/examples/ex27.py

Lines 124 to 137 in 4aa87e1

	def inlet_dofs(self):
	return self.basis['inlet'].find_dofs()['inlet'].all()
	@staticmethod
	def parabolic(x, y):
	"""return the plane Poiseuille parabolic inlet profile"""
	return 4 * y * (1. - y), np.zeros_like(y)
	def make_vector(self):
	"""prepopulate solution vector with Dirichlet conditions"""
	uvp = np.zeros(self.basis['u'].N + self.basis['p'].N)
	I = self.inlet_dofs()
	uvp[I] = L2_projection(self.parabolic, self.basis['inlet'], I)
	return uvp

You don't need the L2_projection if you've only got uniform Dirichlet data.

I'm not sure about get_dofs; I should look into how they work and whether the future is get or find or both.

[bhaveshshrimali]

Thanks so much . I will check it out. Yeah I'm no longer using L2_projection in #455 and also the last updated version in #439

[bhaveshshrimali]

I just took a look at find_dofs and it seems to do a bit more than get_dofs. In fact find_dofs accepts a list of dofnames to skip, which is precisely what is needed here, I think.

scikit-fem/skfem/assembly/basis/basis.py

Lines 91 to 101 in 4aa87e1

	def find_dofs(self,
	facets: Dict[str, ndarray] = None,
	skip: List[str] = None) -> Dict[str, Dofs]:
	"""Return global DOF numbers corresponding to a dictionary of facets.
	Facets can be queried from :class:`~skfem.mesh.Mesh` objects:
	>>> m = MeshTri()
	>>> m.refine()
	>>> m.facets_satisfying(lambda x: x[0] == 0)
	array([1, 5])

So incorporating it, the code looks much cleaner:

dofsold = [
    basis["u"].find_dofs({"left":mesh.facets_satisfying(lambda x: x[0] == 0.)}, skip=["u^2", "u^3"]),
    basis["u"].find_dofs({"bottom":mesh.facets_satisfying(lambda x: x[1] == 0.)}, skip=["u^1", "u^3"]),
    basis["u"].find_dofs({"back":mesh.facets_satisfying(lambda x: x[2] == 0.)}, skip=["u^1", "u^2"]),
    basis["u"].find_dofs({"front":mesh.facets_satisfying(lambda x: x[2] == 1.)}, skip=["u^1", "u^2"])
]

dofs = {}
for dof in dofsold:
    dofs.update(dof)

du[dofs["left"].all()] = 0.
du[dofs["bottom"].all()] = 0.
du[dofs["back"].all()] = 0.
du[dofs["front"].all()] = stretch_

I = basis["u"].complement_dofs(dofs)

Now back to debugging #455

Thanks a lot for the help

[kinnala]

Let me add that ea1e72b introduced these Dofs and DofsView objects, returned by get_dofs and find_dofs. Now they pretty much just reflect the previous API but the plan is that if there are any future requirements or helper methods needed, they can be added directly to DofsView.

Note: There is some vague performance issue I encountered with DofsView for large tetrahedral meshes but didn't debug it any further yet, so there might be some possible improvements there.