d7/d95/assemble__matrix__impl_8h_source.html

 // Copyright (C) 2018-2019 Garth N. Wells

 //

 // This file is part of DOLFINx (https://www.fenicsproject.org)

 //

 // SPDX-License-Identifier:    LGPL-3.0-or-later


 #pragma once


 #include "DofMap.h"

 #include "Form.h"

 #include "utils.h"

 #include <dolfinx/fem/FunctionSpace.h>

 #include <dolfinx/graph/AdjacencyList.h>

 #include <dolfinx/la/utils.h>

 #include <dolfinx/mesh/Geometry.h>

 #include <dolfinx/mesh/Mesh.h>

 #include <dolfinx/mesh/Topology.h>

 #include <functional>

 #include <iterator>

 #include <vector>


 namespace dolfinx::fem::impl

 {


 template <typename T>

 void assemble_matrix(

     const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,

                             const std::int32_t*, const T*)>& mat_set_values,

     const Form<T>& a, const xtl::span<const T>& constants,

     const array2d<T>& coeffs, const std::vector<bool>& bc0,

     const std::vector<bool>& bc1);


 template <typename T>

 void assemble_cells(

     const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,

                             const std::int32_t*, const T*)>& mat_set,

     const mesh::Geometry& geometry,

     const xtl::span<const std::int32_t>& active_cells,

     const std::function<void(const xtl::span<T>&,

                              const xtl::span<const std::uint32_t>&,

                              std::int32_t, int)>& dof_transform,

     const graph::AdjacencyList<std::int32_t>& dofmap0, const int bs0,

     const std::function<void(const xtl::span<T>&,

                              const xtl::span<const std::uint32_t>&,

                              std::int32_t, int)>& dof_transform_to_transpose,

     const graph::AdjacencyList<std::int32_t>& dofmap1, const int bs1,

     const std::vector<bool>& bc0, const std::vector<bool>& bc1,

     const std::function<void(T*, const T*, const T*, const double*, const int*,

                              const std::uint8_t*)>& kernel,

     const array2d<T>& coeffs, const xtl::span<const T>& constants,

     const xtl::span<const std::uint32_t>& cell_info)

 {

   // Prepare cell geometry

   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();


   // FIXME: Add proper interface for num coordinate dofs

   const std::size_t num_dofs_g = x_dofmap.num_links(0);

   const xt::xtensor<double, 2>& x_g = geometry.x();


   // Iterate over active cells

   const int num_dofs0 = dofmap0.links(0).size();

   const int num_dofs1 = dofmap1.links(0).size();

   const int ndim0 = bs0 * num_dofs0;

   const int ndim1 = bs1 * num_dofs1;

   std::vector<T> Ae(ndim0 * ndim1);

   const xtl::span<T> _Ae(Ae);

   std::vector<double> coordinate_dofs(3 * num_dofs_g);


   for (std::int32_t c : active_cells)

   {

     // Get cell coordinates/geometry

     auto x_dofs = x_dofmap.links(c);

     for (std::size_t i = 0; i < x_dofs.size(); ++i)

     {

       std::copy_n(xt::row(x_g, x_dofs[i]).begin(), 3,

                   std::next(coordinate_dofs.begin(), 3 * i));

     }


     // Tabulate tensor

     std::fill(Ae.begin(), Ae.end(), 0);

     kernel(Ae.data(), coeffs.row(c).data(), constants.data(),

            coordinate_dofs.data(), nullptr, nullptr);


     dof_transform(_Ae, cell_info, c, ndim1);

     dof_transform_to_transpose(_Ae, cell_info, c, ndim0);


     // Zero rows/columns for essential bcs

     auto dofs0 = dofmap0.links(c);

     auto dofs1 = dofmap1.links(c);

     if (!bc0.empty())

     {

       for (int i = 0; i < num_dofs0; ++i)

       {

         for (int k = 0; k < bs0; ++k)

         {

           if (bc0[bs0 * dofs0[i] + k])

           {

             // Zero row bs0 * i + k

             const int row = bs0 * i + k;

             std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0.0);

           }

         }

       }

     }


     if (!bc1.empty())

     {

       for (int j = 0; j < num_dofs1; ++j)

       {

         for (int k = 0; k < bs1; ++k)

         {

           if (bc1[bs1 * dofs1[j] + k])

           {

             // Zero column bs1 * j + k

             const int col = bs1 * j + k;

             for (int row = 0; row < ndim0; ++row)

               Ae[row * ndim1 + col] = 0.0;

           }

         }

       }

     }


     mat_set(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(), Ae.data());

   }

 }


 template <typename T>

 void assemble_exterior_facets(

     const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,

                             const std::int32_t*, const T*)>& mat_set,

     const mesh::Mesh& mesh, const xtl::span<const std::int32_t>& active_facets,

     const std::function<void(const xtl::span<T>&,

                              const xtl::span<const std::uint32_t>&,

                              std::int32_t, int)>& dof_transform,

     const graph::AdjacencyList<std::int32_t>& dofmap0, int bs0,

     const std::function<void(const xtl::span<T>&,

                              const xtl::span<const std::uint32_t>&,

                              std::int32_t, int)>& dof_transform_to_transpose,

     const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,

     const std::vector<bool>& bc0, const std::vector<bool>& bc1,

     const std::function<void(T*, const T*, const T*, const double*, const int*,

                              const std::uint8_t*)>& kernel,

     const array2d<T>& coeffs, const xtl::span<const T>& constants,

     const xtl::span<const std::uint32_t>& cell_info,

     const std::function<std::uint8_t(std::size_t)>& get_perm)

 {

   const int tdim = mesh.topology().dim();


   // Prepare cell geometry

   const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();


   // FIXME: Add proper interface for num coordinate dofs

   const std::size_t num_dofs_g = x_dofmap.num_links(0);

   const xt::xtensor<double, 2>& x_g = mesh.geometry().x();


   // Data structures used in assembly

   std::vector<double> coordinate_dofs(3 * num_dofs_g);

   const int num_dofs0 = dofmap0.links(0).size();

   const int num_dofs1 = dofmap1.links(0).size();

   const int ndim0 = bs0 * num_dofs0;

   const int ndim1 = bs1 * num_dofs1;

   std::vector<T> Ae(ndim0 * ndim1);

   const xtl::span<T> _Ae(Ae);


   // Iterate over all facets

   auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);

   assert(f_to_c);

   auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);

   assert(c_to_f);


   for (std::int32_t f : active_facets)

   {

     auto cells = f_to_c->links(f);

     assert(cells.size() == 1);


     // Get local index of facet with respect to the cell

     auto facets = c_to_f->links(cells[0]);

     auto it = std::find(facets.begin(), facets.end(), f);

     assert(it != facets.end());

     const int local_facet = std::distance(facets.begin(), it);


     // Get cell coordinates/geometry

     auto x_dofs = x_dofmap.links(cells[0]);

     for (std::size_t i = 0; i < x_dofs.size(); ++i)

     {

       std::copy_n(xt::row(x_g, x_dofs[i]).begin(), 3,

                   std::next(coordinate_dofs.begin(), 3 * i));

     }


     // Tabulate tensor

     std::uint8_t perm = get_perm(cells[0] * facets.size() + local_facet);

     std::fill(Ae.begin(), Ae.end(), 0);

     kernel(Ae.data(), coeffs.row(cells[0]).data(), constants.data(),

            coordinate_dofs.data(), &local_facet, &perm);


     dof_transform(_Ae, cell_info, cells[0], ndim1);

     dof_transform_to_transpose(_Ae, cell_info, cells[0], ndim0);


     // Zero rows/columns for essential bcs

     auto dofs0 = dofmap0.links(cells[0]);

     auto dofs1 = dofmap1.links(cells[0]);

     if (!bc0.empty())

     {

       for (int i = 0; i < num_dofs0; ++i)

       {

         for (int k = 0; k < bs0; ++k)

         {

           if (bc0[bs0 * dofs0[i] + k])

           {

             // Zero row bs0 * i + k

             const int row = bs0 * i + k;

             std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0.0);

           }

         }

       }

     }

     if (!bc1.empty())

     {

       for (int j = 0; j < num_dofs1; ++j)

       {

         for (int k = 0; k < bs1; ++k)

         {

           if (bc1[bs1 * dofs1[j] + k])

           {

             // Zero column bs1 * j + k

             const int col = bs1 * j + k;

             for (int row = 0; row < ndim0; ++row)

               Ae[row * ndim1 + col] = 0.0;

           }

         }

       }

     }


     mat_set(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(), Ae.data());

   }

 }


 template <typename T>

 void assemble_interior_facets(

     const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,

                             const std::int32_t*, const T*)>& mat_set,

     const mesh::Mesh& mesh, const xtl::span<const std::int32_t>& active_facets,

     const std::function<void(const xtl::span<T>&,

                              const xtl::span<const std::uint32_t>&,

                              std::int32_t, int)>& dof_transform,

     const DofMap& dofmap0, int bs0,

     const std::function<void(const xtl::span<T>&,

                              const xtl::span<const std::uint32_t>&,

                              std::int32_t, int)>& dof_transform_to_transpose,

     const DofMap& dofmap1, int bs1, const std::vector<bool>& bc0,

     const std::vector<bool>& bc1,

     const std::function<void(T*, const T*, const T*, const double*, const int*,

                              const std::uint8_t*)>& kernel,

     const array2d<T>& coeffs, const xtl::span<const int>& offsets,

     const xtl::span<const T>& constants,

     const xtl::span<const std::uint32_t>& cell_info,

     const std::function<std::uint8_t(std::size_t)>& get_perm)

 {

   const int tdim = mesh.topology().dim();


   // Prepare cell geometry

   const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();


   // FIXME: Add proper interface for num coordinate dofs

   const std::size_t num_dofs_g = x_dofmap.num_links(0);

   const xt::xtensor<double, 2>& x_g = mesh.geometry().x();


   // Data structures used in assembly

   xt::xtensor<double, 3> coordinate_dofs({2, num_dofs_g, 3});

   std::vector<T> Ae, be;

   std::vector<T> coeff_array(2 * offsets.back());

   assert(offsets.back() == coeffs.shape[1]);


   // Temporaries for joint dofmaps

   std::vector<std::int32_t> dmapjoint0, dmapjoint1;


   // Iterate over all facets

   auto c = mesh.topology().connectivity(tdim - 1, tdim);

   assert(c);

   auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);

   assert(c_to_f);


   for (std::int32_t facet_index : active_facets)

   {

     // Create attached cells

     auto cells = c->links(facet_index);

     assert(cells.size() == 2);


     // Get local index of facet with respect to the cell

     auto facets0 = c_to_f->links(cells[0]);

     const auto* it0 = std::find(facets0.begin(), facets0.end(), facet_index);

     assert(it0 != facets0.end());

     const int local_facet0 = std::distance(facets0.begin(), it0);

     auto facets1 = c_to_f->links(cells[1]);

     const auto* it1 = std::find(facets1.begin(), facets1.end(), facet_index);

     assert(it1 != facets1.end());

     const int local_facet1 = std::distance(facets1.begin(), it1);


     const std::array local_facet = {local_facet0, local_facet1};


     // Get cell geometry

     auto x_dofs0 = x_dofmap.links(cells[0]);

     for (std::size_t i = 0; i < x_dofs0.size(); ++i)

     {

       std::copy_n(xt::view(x_g, x_dofs0[i]).begin(), 3,

                   xt::view(coordinate_dofs, 0, i, xt::all()).begin());

     }

     auto x_dofs1 = x_dofmap.links(cells[1]);

     for (std::size_t i = 0; i < x_dofs1.size(); ++i)

     {

       std::copy_n(xt::view(x_g, x_dofs1[i]).begin(), 3,

                   xt::view(coordinate_dofs, 1, i, xt::all()).begin());

     }


     // Get dof maps for cells and pack

     xtl::span<const std::int32_t> dmap0_cell0 = dofmap0.cell_dofs(cells[0]);

     xtl::span<const std::int32_t> dmap0_cell1 = dofmap0.cell_dofs(cells[1]);

     dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());

     std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin());

     std::copy(dmap0_cell1.begin(), dmap0_cell1.end(),

               std::next(dmapjoint0.begin(), dmap0_cell0.size()));


     xtl::span<const std::int32_t> dmap1_cell0 = dofmap1.cell_dofs(cells[0]);

     xtl::span<const std::int32_t> dmap1_cell1 = dofmap1.cell_dofs(cells[1]);

     dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());

     std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin());

     std::copy(dmap1_cell1.begin(), dmap1_cell1.end(),

               std::next(dmapjoint1.begin(), dmap1_cell0.size()));


     // Layout for the restricted coefficients is flattened

     // w[coefficient][restriction][dof]

     auto coeff_cell0 = coeffs.row(cells[0]);

     auto coeff_cell1 = coeffs.row(cells[1]);


     // Loop over coefficients

     for (std::size_t i = 0; i < offsets.size() - 1; ++i)

     {

       // Loop over entries for coefficient i

       const int num_entries = offsets[i + 1] - offsets[i];

       std::copy_n(coeff_cell0.data() + offsets[i], num_entries,

                   std::next(coeff_array.begin(), 2 * offsets[i]));

       std::copy_n(coeff_cell1.data() + offsets[i], num_entries,

                   std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));

     }


     const int num_rows = bs0 * dmapjoint0.size();

     const int num_cols = bs1 * dmapjoint1.size();


     // Tabulate tensor

     Ae.resize(num_rows * num_cols);

     std::fill(Ae.begin(), Ae.end(), 0);


     const int facets_per_cell = facets0.size();

     const std::array perm{

         get_perm(cells[0] * facets_per_cell + local_facet[0]),

         get_perm(cells[1] * facets_per_cell + local_facet[1])};

     kernel(Ae.data(), coeff_array.data(), constants.data(),

            coordinate_dofs.data(), local_facet.data(), perm.data());


     dof_transform(Ae, cell_info, cells[0], num_cols);

     dof_transform_to_transpose(Ae, cell_info, cells[0], num_rows);


     // Zero rows/columns for essential bcs

     if (!bc0.empty())

     {

       for (std::size_t i = 0; i < dmapjoint0.size(); ++i)

       {

         for (int k = 0; k < bs0; ++k)

         {

           if (bc0[bs0 * dmapjoint0[i] + k])

           {

             // Zero row bs0 * i + k

             std::fill_n(std::next(Ae.begin(), num_cols * (bs0 * i + k)),

                         num_cols, 0.0);

           }

         }

       }

     }

     if (!bc1.empty())

     {

       for (std::size_t j = 0; j < dmapjoint1.size(); ++j)

       {

         for (int k = 0; k < bs1; ++k)

         {

           if (bc1[bs1 * dmapjoint1[j] + k])

           {

             // Zero column bs1 * j + k

             for (int m = 0; m < num_rows; ++m)

               Ae[m * num_cols + bs1 * j + k] = 0.0;

           }

         }

       }

     }


     mat_set(dmapjoint0.size(), dmapjoint0.data(), dmapjoint1.size(),

             dmapjoint1.data(), Ae.data());

   }

 }


 template <typename T>

 void assemble_matrix(

     const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,

                             const std::int32_t*, const T*)>& mat_set,

     const Form<T>& a, const xtl::span<const T>& constants,

     const array2d<T>& coeffs, const std::vector<bool>& bc0,

     const std::vector<bool>& bc1)

 {

   std::shared_ptr<const mesh::Mesh> mesh = a.mesh();

   assert(mesh);

   const int tdim = mesh->topology().dim();


   // Get dofmap data

   std::shared_ptr<const fem::DofMap> dofmap0

       = a.function_spaces().at(0)->dofmap();

   std::shared_ptr<const fem::DofMap> dofmap1

       = a.function_spaces().at(1)->dofmap();

   assert(dofmap0);

   assert(dofmap1);

   const graph::AdjacencyList<std::int32_t>& dofs0 = dofmap0->list();

   const int bs0 = dofmap0->bs();

   const graph::AdjacencyList<std::int32_t>& dofs1 = dofmap1->list();

   const int bs1 = dofmap1->bs();


   std::shared_ptr<const fem::FiniteElement> element0

       = a.function_spaces().at(0)->element();

   std::shared_ptr<const fem::FiniteElement> element1

       = a.function_spaces().at(1)->element();

   const std::function<void(const xtl::span<T>&,

                            const xtl::span<const std::uint32_t>&, std::int32_t,

                            int)>& dof_transform

       = element0->get_dof_transformation_function<T>();

   const std::function<void(const xtl::span<T>&,

                            const xtl::span<const std::uint32_t>&, std::int32_t,

                            int)>& dof_transform_to_transpose

       = element1->get_dof_transformation_to_transpose_function<T>();


   const bool needs_transformation_data

       = element0->needs_dof_transformations()

         or element1->needs_dof_transformations()

         or a.needs_facet_permutations();

   xtl::span<const std::uint32_t> cell_info;

   if (needs_transformation_data)

   {

     mesh->topology_mutable().create_entity_permutations();

     cell_info = xtl::span(mesh->topology().get_cell_permutation_info());

   }


   for (int i : a.integral_ids(IntegralType::cell))

   {

     const auto& fn = a.kernel(IntegralType::cell, i);

     const std::vector<std::int32_t>& active_cells

         = a.domains(IntegralType::cell, i);

     impl::assemble_cells<T>(mat_set, mesh->geometry(), active_cells,

                             dof_transform, dofs0, bs0,

                             dof_transform_to_transpose, dofs1, bs1, bc0, bc1,

                             fn, coeffs, constants, cell_info);

   }


   if (a.num_integrals(IntegralType::exterior_facet) > 0

       or a.num_integrals(IntegralType::interior_facet) > 0)

   {

     mesh->topology_mutable().create_connectivity(tdim - 1, tdim);

     std::function<std::uint8_t(std::size_t)> get_perm;

     if (a.needs_facet_permutations())

     {

       mesh->topology_mutable().create_entity_permutations();

       const std::vector<std::uint8_t>& perms

           = mesh->topology().get_facet_permutations();

       get_perm = [&perms](std::size_t i) { return perms[i]; };

     }

     else

       get_perm = [](std::size_t) { return 0; };


     for (int i : a.integral_ids(IntegralType::exterior_facet))

     {

       const auto& fn = a.kernel(IntegralType::exterior_facet, i);

       const std::vector<std::int32_t>& active_facets

           = a.domains(IntegralType::exterior_facet, i);

       impl::assemble_exterior_facets<T>(

           mat_set, *mesh, active_facets, dof_transform, dofs0, bs0,

           dof_transform_to_transpose, dofs1, bs1, bc0, bc1, fn, coeffs,

           constants, cell_info, get_perm);

     }


     const std::vector<int> c_offsets = a.coefficient_offsets();

     for (int i : a.integral_ids(IntegralType::interior_facet))

     {

       const auto& fn = a.kernel(IntegralType::interior_facet, i);

       const std::vector<std::int32_t>& active_facets

           = a.domains(IntegralType::interior_facet, i);

       impl::assemble_interior_facets<T>(

           mat_set, *mesh, active_facets, dof_transform, *dofmap0, bs0,

           dof_transform_to_transpose, *dofmap1, bs1, bc0, bc1, fn, coeffs,

           c_offsets, constants, cell_info, get_perm);

     }

   }

 }


 } // namespace dolfinx::fem::impl

dolfinx::array2d
This class provides a dynamic 2-dimensional row-wise array data structure.
Definition: array2d.h:21

dolfinx::array2d::row
constexpr xtl::span< value_type > row(size_type i)
Access a row in the array.
Definition: array2d.h:114

dolfinx::array2d::shape
std::array< size_type, 2 > shape
The shape of the array.
Definition: array2d.h:155

dolfinx::fem::DofMap
Degree-of-freedom map.
Definition: DofMap.h:68

dolfinx::fem::Form
Class for variational forms.
Definition: Form.h:60

dolfinx::fem::Form::integral_ids
std::vector< int > integral_ids(IntegralType type) const
Get the IDs for integrals (kernels) for given integral type. The IDs correspond to the domain IDs whi...
Definition: Form.h:200

dolfinx::fem::Form::kernel
const std::function< void(T *, const T *, const T *, const double *, const int *, const std::uint8_t *)> & kernel(IntegralType type, int i) const
Get the function for 'kernel' for integral i of given type.
Definition: Form.h:161

dolfinx::fem::Form::mesh
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh for the form.
Definition: Form.h:144

dolfinx::fem::Form::needs_facet_permutations
bool needs_facet_permutations() const
Get bool indicating whether permutation data needs to be passed into these integrals.
Definition: Form.h:238

dolfinx::fem::Form::coefficient_offsets
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell. Used to pack data for multiple coefficients in...
Definition: Form.h:243

dolfinx::fem::Form::domains
const std::vector< std::int32_t > & domains(IntegralType type, int i) const
Get the list of mesh entity indices for the ith integral (kernel) for the given domain type,...
Definition: Form.h:217

dolfinx::fem::Form::function_spaces
const std::vector< std::shared_ptr< const fem::FunctionSpace > > & function_spaces() const
Return function spaces for all arguments.
Definition: Form.h:149

dolfinx::fem::Form::num_integrals
int num_integrals(IntegralType type) const
Number of integrals of given type.
Definition: Form.h:186

dolfinx::graph::AdjacencyList
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:47

dolfinx::graph::AdjacencyList::links
xtl::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:130

dolfinx::graph::AdjacencyList::num_links
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:120

dolfinx::mesh::Geometry
Geometry stores the geometry imposed on a mesh.
Definition: Geometry.h:37

dolfinx::mesh::Geometry::dofmap
const graph::AdjacencyList< std::int32_t > & dofmap() const
DOF map.
Definition: Geometry.cpp:22

dolfinx::mesh::Geometry::x
xt::xtensor< double, 2 > & x()
Geometry degrees-of-freedom.
Definition: Geometry.cpp:32

dolfinx::mesh::Mesh
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition: Mesh.h:53

dolfinx::mesh::Mesh::geometry
Geometry & geometry()
Get mesh geometry.
Definition: Mesh.cpp:179

dolfinx::mesh::Mesh::topology
Topology & topology()
Get mesh topology.
Definition: Mesh.cpp:173

dolfinx::mesh::Topology::dim
int dim() const noexcept
Return the topological dimension of the mesh.
Definition: Topology.cpp:163

dolfinx::mesh::Topology::connectivity
std::shared_ptr< const graph::AdjacencyList< std::int32_t > > connectivity(int d0, int d1) const
Return connectivity from entities of dimension d0 to entities of dimension d1.
Definition: Topology.cpp:253

dolfinx::fem::assemble_matrix
void assemble_matrix(const std::function< int(std::int32_t, const std::int32_t *, std::int32_t, const std::int32_t *, const T *)> &mat_add, const Form< T > &a, const xtl::span< const T > &constants, const array2d< T > &coeffs, const std::vector< std::shared_ptr< const DirichletBC< T >>> &bcs)
Assemble bilinear form into a matrix.
Definition: assembler.h:166