Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix extract block for tensor spaces #308

Merged
merged 5 commits into from
Oct 2, 2024
Merged

Conversation

jorgensd
Copy link
Member

@jorgensd jorgensd commented Sep 20, 2024

To be able to easily extract blocks from ufl.MixedFuntionSpaces.

Motivation:
Instead of writing out every block of a and L by hand, as

    a00 = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
    a01 = ufl.inner(c, v) * dx
    a10 = ufl.inner(u, d) * dx

    L0 = ufl.inner(f, v) * dx + ufl.inner(g, v) * ufl.ds
    L1 = ufl.inner(zero, d) * dx

    a = dolfinx.fem.form([[a00, a01], [a10, None]], dtype=dtype)
    L = dolfinx.fem.form([L0, L1], dtype=dtype)

one can now do

    dx = ufl.Measure("dx", domain=mesh)
    f = -ufl.div(ufl.grad(u_ex))
    n = ufl.FacetNormal(mesh)
    g = ufl.dot(ufl.grad(u_ex), n)
    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
    a += ufl.inner(c, v) * dx
    a += ufl.inner(u, d) * dx
    L = ufl.inner(f, v) * dx + ufl.inner(g, v) * ufl.ds
    L += ufl.inner(zero, d) * dx

    a_blocked = ufl.extract_blocks(a)
    L_blocked = ufl.extract_blocks(L)

    a = dolfinx.fem.form(a_blocked, dtype=dtype)
    L = dolfinx.fem.form(L_blocked, dtype=dtype)

It also fixes extract_blocks for tensors, and properly document this function, as to what happens if i or j is not provided.
Now it extracts the blocks as a scalar, vector or tensor.

jorgensd added a commit to scientificcomputing/scifem that referenced this pull request Sep 20, 2024
jorgensd added a commit to scientificcomputing/scifem that referenced this pull request Sep 23, 2024
@@ -148,7 +148,7 @@ def components(self) -> _typing.Dict[_typing.Tuple[int, ...], int]:
offset = 0
c_offset = 0
for e in self.sub_elements:
for i, j in enumerate(np.ndindex(e.value_shape)):
for i, j in enumerate(np.ndindex(e.reference_value_shape)):
Copy link
Member

@mscroggs mscroggs Sep 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This change shouldn't be in this PR (see #307)

Suggested change
for i, j in enumerate(np.ndindex(e.reference_value_shape)):
for i, j in enumerate(np.ndindex(e.value_shape)):

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll revert it once 307 is merged.

@jorgensd jorgensd requested a review from cdaversin September 27, 2024 09:59
Copy link
Contributor

@cdaversin cdaversin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me :)

@jorgensd jorgensd added this pull request to the merge queue Oct 2, 2024
Merged via the queue into main with commit bde4f2f Oct 2, 2024
21 checks passed
@jorgensd jorgensd deleted the dokken/fix-extract-blocks branch October 2, 2024 07:34
@jpdean
Copy link
Member

jpdean commented Oct 2, 2024

Nice! This'll be really handy for coupled problems. Do we have a test?

@jorgensd
Copy link
Member Author

jorgensd commented Oct 2, 2024

Nice! This'll be really handy for coupled problems. Do we have a test?

I think we can adapt the stokes demo.

I use this in many of my own works lately.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants