]> https://gitweb.dealii.org/ - dealii.git/commitdiff
fix a bug
authorJiaqi Zhang <jiaqi2@clemson.edu>
Mon, 15 Aug 2022 05:03:47 +0000 (01:03 -0400)
committerJiaqi Zhang <jiaqi2@clemson.edu>
Wed, 17 Aug 2022 16:09:20 +0000 (12:09 -0400)
include/deal.II/numerics/vector_tools_constraints.templates.h
tests/numerics/no_flux_17.cc
tests/numerics/no_flux_17.output

index 17c05a85eebdfa4850ea9f92ebfc6203cceb94db..14c7c14987cbb82bb7a170bbc5e335fa0a02b5f4 100644 (file)
@@ -461,127 +461,237 @@ namespace VectorTools
             Assert(false, ExcNotImplemented());
         }
     }
-  } // namespace internal
-
-
-  template <int dim, int spacedim>
-  void
-  compute_nonzero_normal_flux_constraints(
-    const DoFHandler<dim, spacedim> &   dof_handler,
-    const unsigned int                  first_vector_component,
-    const std::set<types::boundary_id> &boundary_ids,
-    const std::map<types::boundary_id, const Function<spacedim> *>
-      &                           function_map,
-    AffineConstraints<double> &   constraints,
-    const Mapping<dim, spacedim> &mapping,
-    const IndexSet &              refinement_edge_indices,
-    const unsigned int            level)
-  {
-    Assert(dim > 1,
-           ExcMessage("This function is not useful in 1d because it amounts "
-                      "to imposing Dirichlet values on the vector-valued "
-                      "quantity."));
-
-    const unsigned int mesh_level =
-      (level == numbers::invalid_unsigned_int) ?
-        dof_handler.get_triangulation().n_global_levels() - 1 :
-        level;
 
-    std::vector<types::global_dof_index> face_dofs;
-
-    // create FE and mapping collections for all elements in use by this
-    // DoFHandler
-    const hp::FECollection<dim, spacedim> &fe_collection =
-      dof_handler.get_fe_collection();
-    hp::MappingCollection<dim, spacedim> mapping_collection;
-    for (unsigned int i = 0; i < fe_collection.size(); ++i)
-      mapping_collection.push_back(mapping);
 
-    // TODO: the implementation makes the assumption that all faces have the
-    // same number of dofs
-    AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
-    const unsigned int face_no = 0;
 
-    // now also create a quadrature collection for the faces of a cell. fill
-    // it with a quadrature formula with the support points on faces for each
-    // FE
-    hp::QCollection<dim - 1> face_quadrature_collection;
-    for (unsigned int i = 0; i < fe_collection.size(); ++i)
-      {
-        const std::vector<Point<dim - 1>> &unit_support_points =
-          fe_collection[i].get_unit_face_support_points(face_no);
+    template <int dim, int spacedim>
+    void
+    map_dof_to_normals_on_level(
+      const DoFHandler<dim, spacedim> &   dof_handler,
+      const unsigned int                  first_vector_component,
+      const std::set<types::boundary_id> &boundary_ids,
+      const std::map<types::boundary_id, const Function<spacedim> *>
+        &                              function_map,
+      hp::FEFaceValues<dim, spacedim> &x_fe_face_values,
+      const IndexSet &                 refinement_edge_indices,
+      const unsigned int               level,
+      std::multimap<
+        internal::VectorDoFTuple<dim>,
+        std::pair<Tensor<1, dim>,
+                  typename DoFHandler<dim, spacedim>::cell_iterator>>
+        &dof_to_normals_map,
+      std::map<internal::VectorDoFTuple<dim>, Vector<double>>
+        &dof_vector_to_b_values)
+    {
+      Assert(level < dof_handler.get_triangulation().n_levels(),
+             ExcInternalError());
+
+      std::vector<types::global_dof_index> face_dofs;
+
+      const auto &face_quadrature_collection =
+        x_fe_face_values.get_quadrature_collection();
+
+      // now loop over all cells and all faces
+      std::set<types::boundary_id>::iterator b_id;
+      for (const auto &cell : dof_handler.cell_iterators_on_level(level))
+        if (cell->level_subdomain_id() != numbers::artificial_subdomain_id &&
+            cell->level_subdomain_id() != numbers::invalid_subdomain_id)
+          for (const unsigned int face_no : cell->face_indices())
+            if ((b_id = boundary_ids.find(
+                   cell->face(face_no)->boundary_id())) != boundary_ids.end())
+              {
+                const FiniteElement<dim> &fe = cell->get_fe();
+                typename DoFHandler<dim, spacedim>::level_face_iterator face =
+                  cell->face(face_no);
 
-        Assert(unit_support_points.size() ==
-                 fe_collection[i].n_dofs_per_face(face_no),
-               ExcInternalError());
+                // get the indices of the dofs on this cell...
+                face_dofs.resize(fe.n_dofs_per_face(face_no));
 
-        face_quadrature_collection.push_back(
-          Quadrature<dim - 1>(unit_support_points));
-      }
+                face->get_mg_dof_indices(level,
+                                         face_dofs,
+                                         cell->active_fe_index());
 
-    // now create the object with which we will generate the normal vectors
-    hp::FEFaceValues<dim, spacedim> x_fe_face_values(mapping_collection,
-                                                     fe_collection,
-                                                     face_quadrature_collection,
-                                                     update_quadrature_points |
-                                                       update_normal_vectors);
+                x_fe_face_values.reinit(cell, face_no);
+                const FEFaceValues<dim> &fe_values =
+                  x_fe_face_values.get_present_fe_values();
+
+                // then identify which of them correspond to the selected set of
+                // vector components
+                for (unsigned int i = 0; i < face_dofs.size(); ++i)
+                  if (fe.face_system_to_component_index(i, face_no).first ==
+                      first_vector_component)
+                    // Refinement edge indices are going to be constrained to 0
+                    // during a multigrid cycle and do not need no-normal-flux
+                    // constraints, so skip them:
+                    if (!refinement_edge_indices.is_element(face_dofs[i]))
+                      {
+                        // find corresponding other components of vector
+                        internal::VectorDoFTuple<dim> vector_dofs;
+                        vector_dofs.dof_indices[0] = face_dofs[i];
+
+                        Assert(
+                          first_vector_component + dim <= fe.n_components(),
+                          ExcMessage(
+                            "Error: the finite element does not have enough components "
+                            "to define a normal direction."));
+
+                        for (unsigned int k = 0;
+                             k < fe.n_dofs_per_face(face_no);
+                             ++k)
+                          if ((k != i) &&
+                              (face_quadrature_collection[cell
+                                                            ->active_fe_index()]
+                                 .point(k) ==
+                               face_quadrature_collection[cell
+                                                            ->active_fe_index()]
+                                 .point(i)) &&
+                              (fe.face_system_to_component_index(k, face_no)
+                                 .first >= first_vector_component) &&
+                              (fe.face_system_to_component_index(k, face_no)
+                                 .first < first_vector_component + dim))
+                            vector_dofs.dof_indices
+                              [fe.face_system_to_component_index(k, face_no)
+                                 .first -
+                               first_vector_component] = face_dofs[k];
+
+                        for (unsigned int d = 0; d < dim; ++d)
+                          Assert(vector_dofs.dof_indices[d] <
+                                   dof_handler.n_dofs(),
+                                 ExcInternalError());
+
+                        // we need the normal vector on this face. we know that
+                        // it is a vector of length 1 but at least with higher
+                        // order mappings it isn't always possible to guarantee
+                        // that each component is exact up to zero tolerance. in
+                        // particular, as shown in the deal.II/no_flux_06 test,
+                        // if we just take the normal vector as given by the
+                        // fe_values object, we can get entries in the normal
+                        // vectors of the unit cube that have entries up to
+                        // several times 1e-14.
+                        //
+                        // the problem with this is that this later yields
+                        // constraints that are circular (e.g., in the testcase,
+                        // we get constraints of the form
+                        //
+                        // x22 =  2.93099e-14*x21 + 2.93099e-14*x23
+                        // x21 = -2.93099e-14*x22 + 2.93099e-14*x21
+                        //
+                        // in both of these constraints, the small numbers
+                        // should be zero and the constraints should simply be
+                        // x22 = x21 = 0
+                        //
+                        // to achieve this, we utilize that we know that the
+                        // normal vector has (or should have) length 1 and that
+                        // we can simply set small elements to zero (without
+                        // having to check that they are small *relative to
+                        // something else*). we do this and then normalize the
+                        // length of the vector back to one, just to be on the
+                        // safe side
+                        //
+                        // one more point: we would like to use the "real"
+                        // normal vector here, as provided by the boundary
+                        // description and as opposed to what we get from the
+                        // FEValues object. we do this in the immediately next
+                        // line, but as is obvious, the boundary only has a
+                        // vague idea which side of a cell it is on -- indicated
+                        // by the face number. in other words, it may provide
+                        // the inner or outer normal. by and large, there is no
+                        // harm from this, since the tangential vector we
+                        // compute is still the same. however, we do average
+                        // over normal vectors from adjacent cells and if they
+                        // have recorded normal vectors from the inside once and
+                        // from the outside the other time, then this averaging
+                        // is going to run into trouble. as a consequence we ask
+                        // the mapping after all for its normal vector, but we
+                        // only ask it so that we can possibly correct the sign
+                        // of the normal vector provided by the boundary if they
+                        // should point in different directions. this is the
+                        // case in tests/deal.II/no_flux_11.
+                        Tensor<1, dim> normal_vector =
+                          (cell->face(face_no)->get_manifold().normal_vector(
+                            cell->face(face_no),
+                            fe_values.quadrature_point(i)));
+                        if (normal_vector * fe_values.normal_vector(i) < 0)
+                          normal_vector *= -1;
+                        Assert(std::fabs(normal_vector.norm() - 1) < 1e-14,
+                               ExcInternalError());
+                        for (unsigned int d = 0; d < dim; ++d)
+                          if (std::fabs(normal_vector[d]) < 1e-13)
+                            normal_vector[d] = 0;
+                        normal_vector /= normal_vector.norm();
+
+                        const Point<dim> &point = fe_values.quadrature_point(i);
+                        Vector<double>    b_values(dim);
+                        function_map.at(*b_id)->vector_value(point, b_values);
+
+                        // now enter the (dofs,(normal_vector,cell)) entry into
+                        // the map
+                        dof_to_normals_map.insert(
+                          std::make_pair(vector_dofs,
+                                         std::make_pair(normal_vector, cell)));
+                        dof_vector_to_b_values.insert(
+                          std::make_pair(vector_dofs, b_values));
 
-    // have a map that stores normal vectors for each vector-dof tuple we want
-    // to constrain. since we can get at the same vector dof tuple more than
-    // once (for example if it is located at a vertex that we visit from all
-    // adjacent cells), we will want to average later on the normal vectors
-    // computed on different cells as described in the documentation of this
-    // function. however, we can only average if the contributions came from
-    // different cells, whereas we want to constrain twice or more in case the
-    // contributions came from different faces of the same cell
-    // (i.e. constrain not just the *average normal direction* but *all normal
-    // directions* we find). consequently, we also have to store which cell a
-    // normal vector was computed on
-    using DoFToNormalsMap = std::multimap<
-      internal::VectorDoFTuple<dim>,
-      std::pair<Tensor<1, dim>,
-                typename DoFHandler<dim, spacedim>::cell_iterator>>;
-    std::map<internal::VectorDoFTuple<dim>, Vector<double>>
-      dof_vector_to_b_values;
+#ifdef DEBUG_NO_NORMAL_FLUX
+                        std::cout << "Adding normal vector:" << std::endl
+                                  << "   dofs=" << vector_dofs << std::endl
+                                  << "   cell=" << cell << " at "
+                                  << cell->center() << std::endl
+                                  << "   normal=" << normal_vector << std::endl;
+#endif
+                      }
+              }
+    }
 
-    DoFToNormalsMap dof_to_normals_map;
 
-    // now loop over all cells and all faces
-    std::set<types::boundary_id>::iterator b_id;
-    for (const auto &cell : dof_handler.cell_iterators_on_level(mesh_level))
-      if (cell->level_subdomain_id() != numbers::artificial_subdomain_id &&
-          cell->level_subdomain_id() != numbers::invalid_subdomain_id)
-        for (const unsigned int face_no : cell->face_indices())
-          if ((b_id = boundary_ids.find(cell->face(face_no)->boundary_id())) !=
-              boundary_ids.end())
-            {
-              const FiniteElement<dim> &fe = cell->get_fe();
-              typename DoFHandler<dim, spacedim>::level_face_iterator face =
-                cell->face(face_no);
 
-              // get the indices of the dofs on this cell...
-              face_dofs.resize(fe.n_dofs_per_face(face_no));
+    template <int dim, int spacedim>
+    void
+    map_dof_to_normals(
+      const DoFHandler<dim, spacedim> &   dof_handler,
+      const unsigned int                  first_vector_component,
+      const std::set<types::boundary_id> &boundary_ids,
+      const std::map<types::boundary_id, const Function<spacedim> *>
+        &                              function_map,
+      hp::FEFaceValues<dim, spacedim> &x_fe_face_values,
+      std::multimap<
+        internal::VectorDoFTuple<dim>,
+        std::pair<Tensor<1, dim>,
+                  typename DoFHandler<dim, spacedim>::cell_iterator>>
+        &dof_to_normals_map,
+      std::map<internal::VectorDoFTuple<dim>, Vector<double>>
+        &dof_vector_to_b_values)
+    {
+      std::vector<types::global_dof_index> face_dofs;
+
+      const auto &face_quadrature_collection =
+        x_fe_face_values.get_quadrature_collection();
+
+      // now loop over all cells and all faces
+      std::set<types::boundary_id>::iterator b_id;
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        if (!cell->is_artificial())
+          for (const unsigned int face_no : cell->face_indices())
+            if ((b_id = boundary_ids.find(
+                   cell->face(face_no)->boundary_id())) != boundary_ids.end())
+              {
+                const FiniteElement<dim> &fe = cell->get_fe();
+                typename DoFHandler<dim, spacedim>::face_iterator face =
+                  cell->face(face_no);
 
-              if (level != numbers::invalid_unsigned_int)
-                face->get_mg_dof_indices(mesh_level,
-                                         face_dofs,
-                                         cell->active_fe_index());
-              else
+                // get the indices of the dofs on this cell...
+                face_dofs.resize(fe.n_dofs_per_face(face_no));
                 face->get_dof_indices(face_dofs, cell->active_fe_index());
 
-              x_fe_face_values.reinit(cell, face_no);
-              const FEFaceValues<dim> &fe_values =
-                x_fe_face_values.get_present_fe_values();
+                x_fe_face_values.reinit(cell, face_no);
+                const FEFaceValues<dim> &fe_values =
+                  x_fe_face_values.get_present_fe_values();
 
-              // then identify which of them correspond to the selected set of
-              // vector components
-              for (unsigned int i = 0; i < face_dofs.size(); ++i)
-                if (fe.face_system_to_component_index(i, face_no).first ==
-                    first_vector_component)
-                  // Refinement edge indices are going to be constrained to 0
-                  // during a multigrid cycle and do not need no-normal-flux
-                  // constraints, so skip them:
-                  if (!refinement_edge_indices.is_element(face_dofs[i]))
+                // then identify which of them correspond to the selected set of
+                // vector components
+                for (unsigned int i = 0; i < face_dofs.size(); ++i)
+                  if (fe.face_system_to_component_index(i, face_no).first ==
+                      first_vector_component)
                     {
                       // find corresponding other components of vector
                       internal::VectorDoFTuple<dim> vector_dofs;
@@ -694,7 +804,114 @@ namespace VectorTools
                                 << "   normal=" << normal_vector << std::endl;
 #endif
                     }
-            }
+              }
+    }
+
+  } // namespace internal
+
+
+
+  template <int dim, int spacedim>
+  void
+  compute_nonzero_normal_flux_constraints(
+    const DoFHandler<dim, spacedim> &   dof_handler,
+    const unsigned int                  first_vector_component,
+    const std::set<types::boundary_id> &boundary_ids,
+    const std::map<types::boundary_id, const Function<spacedim> *>
+      &                           function_map,
+    AffineConstraints<double> &   constraints,
+    const Mapping<dim, spacedim> &mapping,
+    const IndexSet &              refinement_edge_indices,
+    const unsigned int            level)
+  {
+    Assert(dim > 1,
+           ExcMessage("This function is not useful in 1d because it amounts "
+                      "to imposing Dirichlet values on the vector-valued "
+                      "quantity."));
+
+    // create FE and mapping collections for all elements in use by this
+    // DoFHandler
+    const hp::FECollection<dim, spacedim> &fe_collection =
+      dof_handler.get_fe_collection();
+    hp::MappingCollection<dim, spacedim> mapping_collection;
+    for (unsigned int i = 0; i < fe_collection.size(); ++i)
+      mapping_collection.push_back(mapping);
+
+    // TODO: the implementation makes the assumption that all faces have the
+    // same number of dofs
+    AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
+    const unsigned int face_no = 0;
+
+    // now also create a quadrature collection for the faces of a cell. fill
+    // it with a quadrature formula with the support points on faces for each
+    // FE
+    hp::QCollection<dim - 1> face_quadrature_collection;
+    for (unsigned int i = 0; i < fe_collection.size(); ++i)
+      {
+        const std::vector<Point<dim - 1>> &unit_support_points =
+          fe_collection[i].get_unit_face_support_points(face_no);
+
+        Assert(unit_support_points.size() ==
+                 fe_collection[i].n_dofs_per_face(face_no),
+               ExcInternalError());
+
+        face_quadrature_collection.push_back(
+          Quadrature<dim - 1>(unit_support_points));
+      }
+
+    // now create the object with which we will generate the normal vectors
+    hp::FEFaceValues<dim, spacedim> x_fe_face_values(mapping_collection,
+                                                     fe_collection,
+                                                     face_quadrature_collection,
+                                                     update_quadrature_points |
+                                                       update_normal_vectors);
+
+    // have a map that stores normal vectors for each vector-dof tuple we want
+    // to constrain. since we can get at the same vector dof tuple more than
+    // once (for example if it is located at a vertex that we visit from all
+    // adjacent cells), we will want to average later on the normal vectors
+    // computed on different cells as described in the documentation of this
+    // function. however, we can only average if the contributions came from
+    // different cells, whereas we want to constrain twice or more in case the
+    // contributions came from different faces of the same cell
+    // (i.e. constrain not just the *average normal direction* but *all normal
+    // directions* we find). consequently, we also have to store which cell a
+    // normal vector was computed on
+    using DoFToNormalsMap = std::multimap<
+      internal::VectorDoFTuple<dim>,
+      std::pair<Tensor<1, dim>,
+                typename DoFHandler<dim, spacedim>::cell_iterator>>;
+    std::map<internal::VectorDoFTuple<dim>, Vector<double>>
+      dof_vector_to_b_values;
+
+    DoFToNormalsMap dof_to_normals_map;
+
+    if (level == numbers::invalid_unsigned_int)
+      {
+        // active cells
+        internal::map_dof_to_normals<dim, spacedim>(dof_handler,
+                                                    first_vector_component,
+                                                    boundary_ids,
+                                                    function_map,
+                                                    x_fe_face_values,
+                                                    dof_to_normals_map,
+                                                    dof_vector_to_b_values);
+      }
+    else
+      { // level cells
+        internal::map_dof_to_normals_on_level<dim, spacedim>(
+          dof_handler,
+          first_vector_component,
+          boundary_ids,
+          function_map,
+          x_fe_face_values,
+          refinement_edge_indices,
+          level,
+          dof_to_normals_map,
+          dof_vector_to_b_values);
+      }
+
+
 
     // Now do something with the collected information. To this end, loop
     // through all sets of pairs (dofs,normal_vector) and identify which
index 9acf3076fad658f38551f131c2e989e2bf812412..edd36d3a766183a801cd96a5219bb5e07980ca50 100644 (file)
@@ -43,7 +43,12 @@ compare_constraints(const AffineConstraints<double> &constraints1,
 {
   const double tol = 1.e-15;
   if (constraints1.n_constraints() != constraints2.n_constraints())
-    return false;
+    {
+      deallog << "n_constraints_1 = " << constraints1.n_constraints()
+              << "n_constraints_2 = " << constraints2.n_constraints()
+              << std::endl;
+      return false;
+    }
 
   auto line1    = (constraints1.get_lines()).begin();
   auto line2    = (constraints2.get_lines()).begin();
@@ -54,7 +59,11 @@ compare_constraints(const AffineConstraints<double> &constraints1,
       const unsigned int line_index  = line1->index;
       const unsigned int line_index2 = line2->index;
       if (line_index != line_index2)
-        return false;
+        {
+          deallog << " line_index_1 = " << line_index
+                  << " line_index_2 = " << line_index2 << std::endl;
+          return false;
+        }
 
       const std::vector<std::pair<types::global_dof_index, double>>
         *constraint_entries_1 = constraints1.get_constraint_entries(line_index);
@@ -132,6 +141,9 @@ get_constraints_on_active_cells(
   VectorTools::compute_no_normal_flux_constraints(
     dofh, 0, no_normal_flux_boundaries, constraints, mapping);
   constraints.close();
+  // deallog<<" ***active cell constraints:
+  // "<<constraints.n_constraints()<<std::endl;
+  // constraints.print(deallog.get_file_stream());
 }
 
 template <int dim>
@@ -169,6 +181,9 @@ get_constraints_on_levels(
                                                       level);
       user_level_constraints.close();
       level_constraints[level].copy_from(user_level_constraints);
+      // deallog<<" ***level cell constraints:
+      // "<<user_level_constraints.n_constraints()<<std::endl;
+      // user_level_constraints.print(deallog.get_file_stream());
     }
 }
 
@@ -306,4 +321,31 @@ main()
         run(triangulation, no_flux_boundary);
       }
   }
+  {
+    // test adaptive refinement
+    // the output file is generated by the old version of
+    // compute_no_normal_flux_constraints() and copied here.
+    deallog << " adaptive hyper_shell:" << std::endl;
+
+    const unsigned int dim = 2;
+    FESystem<dim>      fe(FE_Q<dim>(2), dim);
+
+    const MappingQ<dim> mapping(4);
+    Triangulation<dim>  triangulation(
+      Triangulation<dim>::limit_level_difference_at_vertices);
+    GridGenerator::hyper_shell(triangulation, Point<dim>(), 0.5, 1.);
+    std::set<types::boundary_id> no_flux_boundary{0};
+    triangulation.begin_active()->set_refine_flag();
+    triangulation.execute_coarsening_and_refinement();
+
+    DoFHandler<dim> dofh(triangulation);
+    dofh.distribute_dofs(fe);
+
+    AffineConstraints<double> constraints;
+
+    VectorTools::compute_no_normal_flux_constraints(
+      dofh, 0, no_flux_boundary, constraints, mapping);
+
+    constraints.print(deallog.get_file_stream());
+  }
 }
index 6196d0fe2daafe8708459432a7f9e728628da47a..7d7b0df189c866451c6f038dcc6403e9efc5d43c 100644 (file)
@@ -58,3 +58,48 @@ DEAL::Level 1 OK
 DEAL:: dim 3
 DEAL::Level 0 OK
 DEAL::Level 1 OK
+DEAL:: adaptive hyper_shell:
+    0 1:  -0.7265425
+    3 2:  -0.3249197
+    4 5:  -0.7265425
+    7 6:  -0.3249197
+    13 12:  -0.7265425
+    15 14:  -0.7265425
+    19 18:  0.3249197
+    21 20:  0.3249197
+    25 = 0
+    27 = 0
+    30 31:  0.7265425
+    32 33:  0.7265425
+    37 36:  0.7265425
+    39 38:  0.7265425
+    42 = 0
+    44 = 0
+    48 49:  0.3249197
+    50 51:  0.3249197
+    54 55:  -0.7265425
+    56 57:  -0.7265425
+    60 61:  -0.3249197
+    62 63:  -0.3249197
+    67 66:  -0.3249197
+    69 68:  -0.3249197
+    73 72:  -0.7265425
+    75 74:  -0.7265425
+    79 78:  0.3249197
+    81 80:  0.3249197
+    85 = 0
+    87 = 0
+    90 91:  0.7265425
+    92 93:  0.7265425
+    97 96:  0.7265425
+    99 98:  0.7265425
+    102 = 0
+    104 = 0
+    108 109:  0.3249197
+    110 111:  0.3249197
+    114 115:  -0.3249197
+    124 125:  -0.1583844
+    134 135:  -0.5095254
+    140 141:  -0.3249197
+    146 147:  -0.1583844
+    152 153:  -0.5095254

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.