From 86caace9cab45956557927ac857c1fc2811bc7af Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 18 Jun 2025 15:50:59 +1000 Subject: [PATCH 1/5] Expanded geometry tutorial --- README.md | 16 ++---- src/geometry_dev.jl | 115 +++++++++++++++++++++++++++++++++++++++----- 2 files changed, 108 insertions(+), 23 deletions(-) diff --git a/README.md b/README.md index 300935be..c9fa707f 100644 --- a/README.md +++ b/README.md @@ -2,20 +2,16 @@ *Start solving PDEs in Julia* - | **Documentation** | |:------------ | | [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://gridap.github.io/Tutorials/stable) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://gridap.github.io/Tutorials/dev) | |**Build Status** | | [![Build Status](https://github.com/gridap/Tutorials/workflows/CI/badge.svg?branch=master)](https://github.com/gridap/Tutorials/actions?query=workflow%3ACI) | | **Community** | -| [![Join the chat at https://gitter.im/Gridap-jl/community](https://badges.gitter.im/Gridap-jl/community.svg)](https://gitter.im/Gridap-jl/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) | +| [![Join the chat at https://gitter.im/Gridap-jl/community](https://badges.gitter.im/Gridap-jl/community.svg)](https://gitter.im/Gridap-jl/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) [![Join the chat in Slack](https://img.shields.io/badge/chat-on_slack-purple?style=flat&logo=slack&logoColor=white)(https://julialang.slack.com/archives/C02AK382QSG)]| | **Citation** | | [![DOI](https://joss.theoj.org/papers/10.21105/joss.02520/status.svg)](https://doi.org/10.21105/joss.02520) | - - - ## What This repo contains a set of tutorials to learn how to solve partial differential equations (PDEs) in Julia with the Gridap ecosystem of packages. At the root of this ecosystem is the [Gridap.jl](https://github.com/gridap/Gridap.jl) library. The initial tutorials illustrate the usage of the tools in [Gridap.jl](https://github.com/gridap/Gridap.jl), and these are the tutorials new users should focus on. @@ -35,7 +31,7 @@ Visit one of the following pages, depending of your needs, and start enjoying! - [**STABLE**](https://gridap.github.io/Tutorials/stable) — **Tutorials for the most recently tagged version of Gridap.jl.** - [**DEVEL**](https://gridap.github.io/Tutorials/dev) — *Tutorials for the in-development version of Gridap.jl.* -## Generating tutorials locally +## Generating tutorials locally Note: **only if you intend to contribute to the tutorials as an advanced user/developer** @@ -55,12 +51,11 @@ and then, each time that you perform a change on the tutorial sources, you have julia --project=docs docs/make.jl # From the Unix shell, located at the root of Tutorials repo ``` -to generate the tutorials. The files generated are available at `Tutorials/docs/build/`. - +to generate the tutorials. The files generated are available at `Tutorials/docs/build/`. ## Gridap community -Join to our [gitter](https://gitter.im/Gridap-jl/community) chat to ask questions and interact with the Gridap community. +Join to our [gitter](https://gitter.im/Gridap-jl/community) or [slack](https://julialang.slack.com/archives/C02AK382QSG) chats to ask questions and interact with the Gridap community. ## How to cite Gridap @@ -83,7 +78,4 @@ In order to give credit to the `Gridap` contributors, we simply ask you to cite ## Contact - Please, contact the project administrators, [Santiago Badia](mailto:santiago.badia@monash.edu) and [Francesc Verdugo](mailto:fverdugo@cimne.upc.edu), for further questions about licenses and terms of use. - - diff --git a/src/geometry_dev.jl b/src/geometry_dev.jl index 100fcf80..e356c799 100644 --- a/src/geometry_dev.jl +++ b/src/geometry_dev.jl @@ -26,14 +26,14 @@ using Plots # We begin by defining helper functions that will be essential throughout this tutorial. # These functions help us visualize and work with our mesh structures. -# Convert a CartesianDiscreteModel to an UnstructuredDiscreteModel for more generic handling +# Convert a `CartesianDiscreteModel` to an `UnstructuredDiscreteModel` for more generic handling. function cartesian_model(args...; kwargs...) UnstructuredDiscreteModel(CartesianDiscreteModel(args...; kwargs...)) end -# Visualization function to plot nodes with their IDs -# Input: node_coords - Array of node coordinates -# node_ids - Array of corresponding node IDs +# Visualization function to plot nodes with their IDs. Input: +# - node_coords: Array of node coordinates. +# - node_ids: Array of corresponding node IDs. function plot_node_numbering(node_coords, node_ids) x = map(c -> c[1], node_coords) y = map(c -> c[2], node_coords) @@ -43,8 +43,8 @@ function plot_node_numbering(node_coords, node_ids) vline!(unique(y), linestyle=:dash, color=:grey) end -# Overloaded method to plot node numbering directly from a model -# This function extracts the necessary information from the model and calls the base plotting function +# Overloaded method to plot node numbering directly from a model. +# This function extracts the necessary information from the model and calls the base plotting function. function plot_node_numbering(model) D = num_cell_dims(model) topo = get_grid_topology(model) @@ -83,7 +83,7 @@ end # # ### Key Concept: Nodes vs. Vertices # -# One of the most important distinctions in Gridap is between nodes and vertices: +# A very important distinction in Gridap is between nodes and vertices: # # - **Vertices** (Topological entities): # * 0-dimensional entities in the `GridTopology` @@ -91,7 +91,7 @@ end # * Used for neighbor queries and mesh traversal # * Number of vertices depends only on topology # -# - **Nodes** (Geometric entities): +# - **Nodes** (Geometrical entities): # * Control points stored in the `Grid` # * Define the geometry of elements # * Used for interpolation and mapping @@ -180,10 +180,10 @@ node_coordinates = get_node_coordinates(grid) # Physical coordinates of nodes # # There are two ways to get the coordinates of nodes for each cell: # -# 1. Using standard Julia mapping: +# A) Using standard Julia mapping: cell_to_node_coords = map(nodes -> node_coordinates[nodes], cell_to_nodes) -# 2. Using Gridap's lazy evaluation system (more efficient for large meshes): +# B) Using Gridap's lazy evaluation system (more efficient for large meshes): cell_to_node_coords = lazy_map(Broadcasting(Reindex(node_coordinates)),cell_to_nodes) # ### Geometric Mappings @@ -230,7 +230,7 @@ writevtk(new_grid,"half_cylinder_linear") # our half-cylinder looks faceted. This is because we're still using linear elements # (straight edges) to approximate the curved geometry. # -# ### Solution: High-order Elements +# ### Example: High-order Elements # # To accurately represent curved geometries, we need high-order elements: @@ -389,3 +389,96 @@ mobius = UnstructuredDiscreteModel(grid,topo,labels) # Visualize the vertex numbering: plot_node_numbering(mobius) # ![](../assets/geometry/mobius.png) + +# ## FaceLabelings and boundary conditions +# +# The `FaceLabeling` component of a `DiscreteModel` is the way Gridap handles boundary conditions. +# The basic idea is that, similar to Gmsh, we classify the d-faces (cells, faces, edges, nodes) of the mesh +# into different entities (physical groups in Gmsh terminology) which in turn have one or more +# tags/labels associated with them. +# We can then query the `FaceLabeling` for the tags associated with a given d-face, +# or the d-faces associated with a given tag. +# +# We will now explore ways to create and manipulate `Facelabeling` objects. +# +# ### Creating FaceLabelings +# +# The simplest way to create a blank `FaceLabeling` is to use your `GridTopology`: +# + +model = cartesian_model((0,1,0,1),(3,3)) +topo = get_grid_topology(model) + +labels = FaceLabeling(topo) + +# The above `FaceLabeling` is by default created with 2 entities and 2 tags, associated to +# interior and boundary d-faces respectively. The boundary facets are chosen as the ones +# with a single neighboring cell. +# +# We can extract the low-level information from the `FaceLabeling` object: + +tag_names = get_tag_name(labels) # Each name is a string +tag_entities = get_tag_entities(labels) # For each tag, a vector of entities +cell_to_entity = get_face_entity(labels,2) # For each cell, its associated entity +edge_to_entity = get_face_entity(labels,1) # For each edge, its associated entity +node_to_entity = get_face_entity(labels,0) # For each node, its associated entity + +# It is usually more convenient to visualise it in Paraview by exporting to vtk: + +writevtk(model,"labels_basic",labels=labels) + +# Another useful way to create a `FaceLabeling` is by providing a coloring for the mesh cells, +# where each color corresponds to a different tag. +# The d-faces of the mesh will have all the tags associated to the cells that share them. + +cell_to_tag = [1,1,1,2,2,3,2,2,3] +tag_to_name = ["A","B","C"] +labels_cw = Geometry.face_labeling_from_cell_tags(topo,cell_to_tag,tag_to_name) +writevtk(model,"labels_cellwise",labels=labels_cw) + +# We can also create a `FaceLabeling` from a vertex filter. The resulting `FaceLabeling` will have +# only one tag, gathering the d-faces whose vertices ALL fullfill `filter(x) == true`. + +filter(x) = abs(x[1]- 1.0) < 1.e-5 +labels_vf = Geometry.face_labeling_from_vertex_filter(topo, "top", filter) +writevtk(model,"labels_filter",labels=labels_vf) + +# `FaceLabeling` objects can also be merged together. The resulting `FaceLabeling` will have +# the union of the tags and entities of the original ones. +# Note that this modifies the first `FaceLabeling` in place. + +labels = merge!(labels, labels_cw, labels_vf) +writevtk(model,"labels_merged",labels=labels) + +# ### Creating new tags from existing ones +# +# Tags in a `FaceLabeling` support all the usual set operation, i.e union, intersection, +# difference and complementary. + +cell_to_tag = [1,1,1,2,2,3,2,2,3] +tag_to_name = ["A","B","C"] +labels = Geometry.face_labeling_from_cell_tags(topo,cell_to_tag,tag_to_name) + +# Union: Takes as input a list of tags and creates a new tag that is the union of all of them. +Geometry.add_tag_from_tags!(labels,"A∪B",["A","B"]) + +# Intersection: Takes as input a list of tags and creates a new tag that is the intersection of all of them. +Geometry.add_tag_from_tags_intersection!(labels,"A∩B",["A","B"]) + +# Complementary: Takes as input a list of tags and creates a new tag that is the complementary of the union. +Geometry.add_tag_from_tags_complementary!(labels,"!A",["A"]) + +# Set difference: Takes as input two lists of tags (tags_include - tags_exclude) +# and creates a new tag that contains all the d-faces that are in the first list but not in the second. +Geometry.add_tag_from_tags_setdiff!(labels,"A-B",["A"],["B"]) # set difference + +writevtk(model,"labels_setops",labels=labels) + +# ### FaceLabeling queries +# +# The most common way of query information from a `FaceLabeling` is to query a face mask for +# a given tag and face dimension. If multiple tags are provided, the union of the tags is returned. + +face_dim = 1 +mask = get_face_mask(labels,["A","C"],face_dim) # Boolean mask +ids = findall(edges_mask_A) # Edge IDs From 658fbb82e4a08d598bcd2f89bdb94305fb2ce9aa Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 18 Jun 2025 16:10:09 +1000 Subject: [PATCH 2/5] Minor --- src/geometry_dev.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/geometry_dev.jl b/src/geometry_dev.jl index e356c799..7a3fcd06 100644 --- a/src/geometry_dev.jl +++ b/src/geometry_dev.jl @@ -6,6 +6,7 @@ # - How to extract geometrical information from a `Grid`. # - How periodicity is handled in Gridap, and the difference between nodes and vertices. # - How to create a periodic model from scratch, use the example of a Mobius strip. +# - How to create and manipulate `FaceLabeling` objects, which are used to handle boundary conditions. # # ## Required Packages @@ -21,6 +22,7 @@ using Plots # 4. Geometric Mappings # 5. High-order Grids # 6. Periodicity in Gridap +# 7. FaceLabelings # # ## 1. Utility Functions # We begin by defining helper functions that will be essential throughout this tutorial. @@ -390,7 +392,7 @@ mobius = UnstructuredDiscreteModel(grid,topo,labels) plot_node_numbering(mobius) # ![](../assets/geometry/mobius.png) -# ## FaceLabelings and boundary conditions +# ## 7. FaceLabelings and boundary conditions # # The `FaceLabeling` component of a `DiscreteModel` is the way Gridap handles boundary conditions. # The basic idea is that, similar to Gmsh, we classify the d-faces (cells, faces, edges, nodes) of the mesh From 0cb0340218d9ef64fba1313473b8fdcaefdef166 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 18 Jun 2025 16:16:38 +1000 Subject: [PATCH 3/5] Minor --- src/stokes_blocks.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/stokes_blocks.jl b/src/stokes_blocks.jl index 9327e570..96c7270d 100644 --- a/src/stokes_blocks.jl +++ b/src/stokes_blocks.jl @@ -1,4 +1,3 @@ -# # Incompressible Stokes equations in a 2D/3D cavity # # This example solves the incompressible Stokes equations, given by # From a598923571f0ef18f6636ec19f3e15d250cd8184 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 18 Jun 2025 16:27:00 +1000 Subject: [PATCH 4/5] Minor --- src/geometry_dev.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/geometry_dev.jl b/src/geometry_dev.jl index 7a3fcd06..646602cd 100644 --- a/src/geometry_dev.jl +++ b/src/geometry_dev.jl @@ -441,8 +441,8 @@ writevtk(model,"labels_cellwise",labels=labels_cw) # We can also create a `FaceLabeling` from a vertex filter. The resulting `FaceLabeling` will have # only one tag, gathering the d-faces whose vertices ALL fullfill `filter(x) == true`. -filter(x) = abs(x[1]- 1.0) < 1.e-5 -labels_vf = Geometry.face_labeling_from_vertex_filter(topo, "top", filter) +vfilter(x) = abs(x[1]- 1.0) < 1.e-5 +labels_vf = Geometry.face_labeling_from_vertex_filter(topo, "top", vfilter) writevtk(model,"labels_filter",labels=labels_vf) # `FaceLabeling` objects can also be merged together. The resulting `FaceLabeling` will have From 7e6f4476b1683e1e771031e8ce32eec3f32f4096 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 18 Jun 2025 16:39:50 +1000 Subject: [PATCH 5/5] Minor --- src/geometry_dev.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geometry_dev.jl b/src/geometry_dev.jl index 646602cd..68f7c7b4 100644 --- a/src/geometry_dev.jl +++ b/src/geometry_dev.jl @@ -483,4 +483,4 @@ writevtk(model,"labels_setops",labels=labels) face_dim = 1 mask = get_face_mask(labels,["A","C"],face_dim) # Boolean mask -ids = findall(edges_mask_A) # Edge IDs +ids = findall(mask) # Edge IDs