test Namespace Reference

Namespace to group the methods used in testing of the library. More...

Data Structures
struct	testNSearchData

Functions
double	getExactIntegrationRefTri (size_t alpha, size_t beta)
	Computes integration of polynomial exactly over reference triangle.
double	getExactIntegrationRefQuad (size_t alpha, size_t beta)
	Computes integration of polynomial exactly over reference quadrangle.
double	getExactIntegrationRefTet (size_t alpha, size_t beta, size_t theta)
	Computes integration of polynomial exactly over reference tetrahedral.
double	getNChooseR (size_t n, size_t r)
	Computes \( {n\choose r}\) "n choose r".
void	testGraphPartitioningSimple ()
	Tests metis partitioning of graph.
void	testGraphPartitioning (size_t nPart=4, size_t nGrid=10, size_t mHorizon=3, size_t testOption=0, std::string meshFilename="")
	Tests metis partitioning of graph from a 2-D mesh with nonlocal interaction.
template<int dim = 3>
std::string	testNanoflann (size_t N, double L, double dL, int seed)
	Perform test on nsearch.
template<int dim = 3>
std::string	testNanoflannExcludeInclude (size_t N, double L, double dL, int seed, testNSearchData &data)
	Perform test on nsearch.
std::string	testNanoflannClosestPoint (size_t N, double L, double dL, int seed)
	Perform test on nsearch.
std::string	testTaskflow (size_t N, int seed)
	Perform test on taskflow.
void	testMPI (size_t nGrid=10, size_t mHorizon=3, size_t testOption=0, std::string meshFilename="")
	Perform parallelization test using MPI on mesh partition based on metis.
void	testUtilMethods ()
	Test methods
void	testContactStiffness ()
	Checks the contact-stiffness helpers are bit-for-bit what the expressions they replaced produced.
void	testLineElem (size_t n, std::string filepath)
	Perform test on quadrature points on line elements (NOT IMPLEMENTED)
void	testTriElem (size_t n, std::string filepath)
	Perform test on quadrature points on triangle elements.
void	testQuadElem (size_t n, std::string filepath)
	Perform test on quadrature points on quadrangle elements.
void	testTetElem (size_t n, std::string filepath)
	Perform test on quadrature points on tetrahedral elements.
void	testTriElemTime (size_t n, size_t N)
	Computes the time needed when quad data for elements are stored and when they are computed as and when needed.
void	testPatchTri ()
	Perform test on quadrature points on line elements (NOT IMPLEMENTED)
void	testPatchQuad ()
	Perform test on quadrature points on line elements (NOT IMPLEMENTED)
void	testPatchTet ()
	Perform test on quadrature points on line elements (NOT IMPLEMENTED)
void	testPatchTriDistorted ()
	Perform test on quadrature points on line elements (NOT IMPLEMENTED)

Detailed Description

Namespace to group the methods used in testing of the library.

Function Documentation

getExactIntegrationRefQuad()

double test::getExactIntegrationRefQuad	(	size_t	alpha,
		size_t	beta
	)

Computes integration of polynomial exactly over reference quadrangle.

Given \( f(s,t) = s^\alpha\, t^\beta \), the exact integration is given by

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds. \]

If either \( \alpha\) or \( \beta\) are odd number then \( I_{exact} = 0\). Otherwise, \( I_{exact} = \frac{4}{(\alpha +1) (\beta+1)} \).

Parameters

alpha	Polynomial order in variable s
beta	Polynomial order in variable t

Returns

I Exact integration of \( f(s,t) = s^\alpha\, t^\beta \)

Definition at line 194 of file testFeLib.cpp.

                                                                 {
 
  // compute exact integration of s^\alpha t^\beta
  if (alpha % 2 == 0 and beta % 2 == 0)
    return 4. / double((alpha + 1) * (beta + 1));
  else
    return 0.;
}

Referenced by testQuadElem().

Here is the caller graph for this function:

getExactIntegrationRefTet()

double test::getExactIntegrationRefTet	(	size_t	alpha,
		size_t	beta,
		size_t	theta
	)

Computes integration of polynomial exactly over reference tetrahedral.

Parameters

alpha	Polynomial order in variable s
beta	Polynomial order in variable t
theta	Polynomial order in variable r

Returns

I Exact integration of \( f(s,t) = s^\alpha\, t^\beta \, r^\theta \)

Definition at line 203 of file testFeLib.cpp.

                                                     {
 
  double I = 0.;
  for (size_t i = 0; i <= theta + 1; i++) {
 
    double factor_i = test::getNChooseR(theta + 1, i) /
                      (double(theta + 1) * double(i + beta + 1));
    if (i % 2 != 0)
      factor_i = factor_i * (-1.);
 
    for (size_t j = 0; j <= theta + beta + 2 + 1; j++) {
 
      double factor_j =
          test::getNChooseR(theta + beta + 2, j) / (double(j + alpha + 1));
      if (j % 2 != 0)
        factor_j = factor_j * (-1.);
 
      I += factor_i * factor_j;
    }
  }
 
  return I;
}

References getNChooseR().

Referenced by testTetElem().

Here is the call graph for this function:

Here is the caller graph for this function:

getExactIntegrationRefTri()

double test::getExactIntegrationRefTri	(	size_t	alpha,
		size_t	beta
	)

Computes integration of polynomial exactly over reference triangle.

Given \( f(s,t) = s^\alpha\, t^\beta \), the exact integration is given by

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds = \sum_{i=0}^{\beta+1} (-1)^i \frac{{{\beta + 1} \choose i}}{(\alpha + i +1) (\beta + 1)}, \]

where

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

We have \( {a \choose 0} = 1 \) so that term for \( i=0\) is not zero. Above formula gives the exact value of integral of \( f(s,t) = s^\alpha\, t^\beta \) over reference triangle.

Parameters

alpha	Polynomial order in variable s
beta	Polynomial order in variable t

Returns

I Exact integration of \( f(s,t) = s^\alpha\, t^\beta \)

Definition at line 178 of file testFeLib.cpp.

                                                                {
 
  // compute exact integration of s^\alpha t^\beta
  double I = 0.;
  for (size_t k = 0; k <= beta + 1; k++) {
    if (k % 2 == 0)
      I +=
          test::getNChooseR(beta + 1, k) / double((alpha + 1 + k) * (beta + 1));
    else
      I -=
          test::getNChooseR(beta + 1, k) / double((alpha + 1 + k) * (beta + 1));
  }
 
  return I;
}

References getNChooseR().

Referenced by testTriElem().

Here is the call graph for this function:

Here is the caller graph for this function:

getNChooseR()

double test::getNChooseR	(	size_t	n,
		size_t	r
	)

Computes \( {n\choose r}\) "n choose r".

Computes formula

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

Parameters

n	Number
r	Number which is smaller or equal to n

Returns

Value Value of "n choose r"

Definition at line 166 of file testFeLib.cpp.

                                           {
 
  if (r == 0)
    return 1.;
 
  double a = 1.;
  for (size_t i = 1; i <= r; i++)
    a *= double(n - i + 1) / double(i);
 
  return a;
}

Referenced by getExactIntegrationRefTet(), and getExactIntegrationRefTri().

Here is the caller graph for this function:

testContactStiffness()

void test::testContactStiffness

(

)

Checks the contact-stiffness helpers are bit-for-bit what the expressions they replaced produced.

util::normalContactStiffness and util::selfContactStiffness replaced the same formula at nine call sites. The two differ only in the order of operations, which changes the result in the last bit, so each has to reproduce the grouping of the expression it replaced. Otherwise collecting them changes the contact force in every example.

Definition at line 115 of file testUtilLib.cpp.

                                {
 
  // Bulk moduli and horizons that appear in the examples and tests.
  const std::vector<double> Ks = {2.16e7, 1.0e4,     1.0e5,
                                  159.2e9, 216000.0, 2.0e9, 1.23e9};
  const std::vector<double> hs = {6.0e-4, 2.0e-4, 4.0e-4,
                                  3.0e-3, 3.2e-4, 1.25e-3, 3.75e-3};
 
  size_t n_checked = 0;
  for (auto K1 : Ks) {
    for (auto K2 : Ks) {
      for (auto h : hs) {
        // The form the example and test drivers used, at both exponents.
        for (int p : {4, 5}) {
          const double expected =
              18.0 * util::harmonicMean(K1, K2) / (M_PI * std::pow(h, p));
          const double got = util::normalContactStiffness(K1, K2, h, p);
          if (got != expected)
            errExit(std::format(
                "normalContactStiffness({}, {}, {}, {}) = {:.20g}, expected "
                "{:.20g}: the replaced expression gave a different value\n",
                K1, K2, h, p, got, expected));
          ++n_checked;
        }
      }
    }
    for (auto h : hs) {
      // The form BaseParticle uses for internal contact. The grouping is
      // (18 / (pi h^5)) * K and not 18 K / (pi h^5).
      const double expected = (18. / (M_PI * std::pow(h, 5))) * K1;
      const double got = util::selfContactStiffness(K1, h);
      if (got != expected)
        errExit(std::format(
            "selfContactStiffness({}, {}) = {:.20g}, expected {:.20g}: the "
            "order of operations differs from the one in BaseParticle\n",
            K1, h, got, expected));
      ++n_checked;
    }
  }
 
  // The two forms are the same quantity because harmonicMean(K, K) == K.
  for (auto K : Ks)
    if (util::harmonicMean(K, K) != K)
      errExit(std::format("harmonicMean({0}, {0}) != {0}\n", K));
 
  // A horizon of zero divides by zero, so it is rejected instead.
  bool threw = false;
  try {
    util::normalContactStiffness(1., 1., 0.);
  } catch (const std::exception &) {
    threw = true;
  }
  if (!threw)
    errExit("normalContactStiffness accepted a zero horizon\n");
 
  std::cout << std::format(
      "testContactStiffness: {} values equal to the replaced expressions\n",
      n_checked);
}

References anonymous_namespace{testUtilLib.cpp}::errExit(), util::harmonicMean(), util::normalContactStiffness(), and util::selfContactStiffness().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testGraphPartitioning()

void test::testGraphPartitioning	(	size_t	nPart = 4,
		size_t	nGrid = 10,
		size_t	mHorizon = 3,
		size_t	testOption = 0,
		std::string	meshFilename = ""
	)

Tests metis partitioning of graph from a 2-D mesh with nonlocal interaction.

Parameters

nPart	Number of partitions
nGrid	Number of element along a line (total number of elements is N*N)
mHorizon	Integer factor that is used to compute nonlocal radius, i.e., horizon (epsilon = m * h, h being mesh size)
testOption	Test otion flag. 0 - use in-built uniform mesh, 1 - use user-specified mesh
meshFilename	Mesh filename with relative path from the directory where test command is run

Definition at line 137 of file testMeshPartitioningLib.cpp.

                                                                                                                       {
 
  std::cout << "\nMETIS_TEST\n";
  std::cout << "\n  Test the METIS library for graph partitioning for realistic mesh with nonlocal interaction.\n" ;
  std::cout << std::format("\n  Arguments: nPart = {}, nGrid = {}, mHorizon = {}\n", nPart, nGrid, mHorizon);
 
  // create uniform mesh on domain [0, Lx]x[0, Ly]
  auto t1 = steady_clock::now();
 
  // empty mesh object
  size_t dim(2);
  auto mesh = mesh::Mesh(dim);
  mesh.d_spatialDiscretization = "finite_difference";
 
  // create mesh
  std::string outMeshFilename = "";
  if (testOption == 1) {
    // set geometry details
    std::pair<std::vector<double>, std::vector<double>> box;
    std::vector<size_t> nGridVec;
    for (size_t i=0; i<dim; i++) {
      box.first.push_back(0.);
      box.second.push_back(1.);
      nGridVec.push_back(nGrid);
    }
 
    // call utility function to create mesh
    mesh::createUniformMesh(&mesh, dim, box, nGridVec);
 
    // filename for outputting
    outMeshFilename = std::format("uniform_mesh_Lx_{}_Ly_{}_Nx_{}_Ny_{}",
                                               box.second[0], box.second[1],
                                               nGridVec[0], nGridVec[1]);
  }
  else if (testOption == 2) {
    if (meshFilename.empty()) {
      std::cerr << "testGraphPartitioning(): mesh filename is empty.\n";
      exit(1);
    }
 
    // call in-built function of mesh to create data from file
    mesh.createData(meshFilename);
 
    // find the name of mesh file without path and without extension (i.e., remove .vtu/.msh/.csv extension)
    auto f1 = util::io::getFilenameFromPath(meshFilename);
    outMeshFilename = util::io::removeExtensionFromFile(f1);
  }
  else {
    std::cerr << "testGraphPartitioning() accepts either 0 or 1 for testOption. The value "
              << testOption << " is invalid.\n";
    exit(1);
  }
 
  // set nonlocal lengthscale
  double horizon = mHorizon*mesh.d_h;
 
  // print mesh data and write mesh to a file
  std::ostringstream msg;
  msg << mesh.printStr();
  std::cout << msg.str();
 
  auto t2 = steady_clock::now();
  auto setup_time = util::methods::timeDiff(t1, t2, "microseconds");
  std::cout << std::format("Setup time (ms) = {}. \n", setup_time);
 
  // create neighborhood of each node (to be used in metis partitioning of the graph)
  std::vector<std::vector<size_t>> nodeNeighs(mesh.d_numNodes);
  geom::computeNonlocalNeighborhood(mesh.d_nodes, horizon, nodeNeighs);
  auto t3 = steady_clock::now();
  auto neigh_time = util::methods::timeDiff(t2, t3, "microseconds");
  std::cout << std::format("Neighborhood calculation time (ms) = {}.\n", neigh_time);
 
  // at this stage, we have mesh and nonlocal neighborhood
  // we are ready to cast the nonlocal neighborhood into graph and call metis for partitioning of nodes
  std::vector<size_t> nodePartitionRecursive(mesh.d_numNodes, 0);
  std::vector<size_t> nodePartitionKWay(mesh.d_numNodes, 0);
 
  // recursive method
  auto t4 = steady_clock::now();
  mesh::metisGraphPartition("metis_recursive", nodeNeighs, nodePartitionRecursive, nPart);
  auto t5 = steady_clock::now();
 
  // K-way method
  mesh::metisGraphPartition("metis_kway", nodeNeighs, nodePartitionKWay, nPart);
  auto t6 = steady_clock::now();
 
  auto partition_recursive_time = util::methods::timeDiff(t4, t5, "microseconds");
  auto partition_kway_time = util::methods::timeDiff(t5, t6, "microseconds");
  std::cout << std::format("Partition (Recursive) calculation time (ms) = {}.\n", partition_recursive_time);
  std::cout << std::format("Partition (KWay) calculation time (ms) = {}.\n", partition_kway_time);
 
  // write data to file
  outMeshFilename = outMeshFilename + std::format("_mHorizon_{}_nPart_{}", mHorizon, nPart);
  std::cout << "out mesh filename = " << outMeshFilename << std::endl;
  auto writer = rw::writer::Writer(outMeshFilename, "vtu");
  writer.appendMesh(&mesh.d_nodes, mesh.d_eType, &mesh.d_enc);
  writer.appendPointData("Nodal_Volume", &mesh.d_vol);
  writer.appendPointData("Nodal_Partition_Metis_Recursive_Index", &nodePartitionRecursive);
  writer.appendPointData("Nodal_Partition_Metis_KWay_Index", &nodePartitionKWay);
  writer.close();
}

References geom::computeNonlocalNeighborhood(), mesh::createUniformMesh(), util::io::getFilenameFromPath(), mesh::metisGraphPartition(), util::io::removeExtensionFromFile(), and util::methods::timeDiff().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testGraphPartitioningSimple()

void test::testGraphPartitioningSimple

(

)

Tests metis partitioning of graph.

Definition at line 126 of file testMeshPartitioningLib.cpp.

                                       {
 
  printf("\n");
  printf("METIS_TEST\n");
  printf("  Test the METIS library for graph partitioning (simple).\n");
 
  // check that metis is linked and returns a partition
  partGraphKwayTestSimple();
  partGraphRecursiveTestSimple();
}

References anonymous_namespace{testMeshPartitioningLib.cpp}::partGraphKwayTestSimple(), and anonymous_namespace{testMeshPartitioningLib.cpp}::partGraphRecursiveTestSimple().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testLineElem()

void test::testLineElem	(	size_t	n,
		std::string	filepath
	)

Perform test on quadrature points on line elements (NOT IMPLEMENTED)

This function performs accuracy test of the quadrature points for integration over reference line with vertices at {-1, 1}. List of tests are as follows:

1. Computes quadrature points of the given order, writes them to the file, and checks if the sum of quadrature weights is equal to 0.5 (area of reference triangle).

Also tests the exactness of the integration of the polynomial upto given order. Suppose \( n\) is the order of quadrature point, then we test if the integration of the function \( f(s,t) = s^\alpha\, t^\beta \) is exact for \( \alpha \) and \( \beta \) such that \( \alpha+\beta \leq n \). The exact integration of function \( f\) over reference triangle is

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds = \sum_{i=0}^{\beta+1} (-1)^i \frac{{{\beta + 1} \choose i}}{(\alpha + i +1) (\beta + 1)}, \]

where

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

We have \( {a \choose 0} = 1 \) so that term for \( i=0\) is not zero. Above formula gives the exact value of integral of \( f(s,t) = s^\alpha\, t^\beta \) over reference triangle. Approximation by quadrature point is as follows

\[ I_{approx} = \sum_{q=1}^{Q} w_q f(s_q, t_q) \]

where \(Q\) is the total number of quad points, \( w_q\) and \((s_q, t_q)\) are the \( q^{th} \) quad weight and point. In this test, we compare \( I_{exact} \) and \( I_{approx} \) and report problem if both do not match.

1. Test the accuracy for simple quadrangle mesh on square domain [0,1]^2. Exact integration of polynomial \( f(x,y) = x^\alpha \, y^\beta \) on \([0,1]^2\) is given by

   \[ I_{exact} = \frac{1}{(\alpha+1) (\beta+1)}. \]

   We compare above with the approximation computed from the quadrature points. We consider \(\alpha + \beta \leq n \) where \( n\) is the order of approximation we are testing.

Parameters

n	Order of quadrature point approximation
filepath	Path where mesh data for test can be found (expects files 'triMesh_nodes.csv' and 'triMesh_elements.csv' inside the filepath)

Definition at line 229 of file testFeLib.cpp.

{ return; }

Referenced by main().

Here is the caller graph for this function:

testMPI()

void test::testMPI	(	size_t	nGrid = 10,
		size_t	mHorizon = 3,
		size_t	testOption = 0,
		std::string	meshFilename = ""
	)

Perform parallelization test using MPI on mesh partition based on metis.

Parameters

nGrid	Number of element along a line (total number of elements is N*N)
mHorizon	Integer factor that is used to compute nonlocal radius, i.e., horizon (epsilon = m * h, h being mesh size)
testOption	Test otion flag. 0 - use in-built uniform mesh, 1 - use user-specified mesh
meshFilename	Mesh filename with relative path from the directory where test command is run

Definition at line 284 of file testParallelCompLib.cpp.

                                                              {
  int mpiSize, mpiRank;
  MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
  MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
 
  // number of partitions
  size_t nPart(mpiSize);
 
  // create uniform mesh
  size_t dim(2);
  auto mesh = mesh::Mesh(dim);
  mesh.d_spatialDiscretization = "finite_difference";
 
  // create mesh
  std::string outMeshFilename = "";
  if (testOption == 1) {
    // set geometry details
    std::pair<std::vector<double>, std::vector<double>> box;
    std::vector<size_t> nGridVec;
    for (size_t i=0; i<dim; i++) {
      box.first.push_back(0.);
      box.second.push_back(1.);
      nGridVec.push_back(nGrid);
    }
 
    // call utility function to create mesh
    util::io::print("\n\nCreating uniform mesh\n\n");
    mesh::createUniformMesh(&mesh, dim, box, nGridVec);
 
    // filename for outputting
    outMeshFilename = std::format("uniform_mesh_Lx_{}_Ly_{}_Nx_{}_Ny_{}",
                                  box.second[0], box.second[1],
                                  nGridVec[0], nGridVec[1]);
  }
  else if (testOption == 2) {
    if (meshFilename.empty()) {
      std::cerr << "testGraphPartitioning(): mesh filename is empty.\n";
      exit(1);
    }
 
    // call in-built function of mesh to create data from file
    util::io::print("\n\nReading mesh\n\n");
    mesh.createData(meshFilename);
 
    // find the name of mesh file excluding path and extension
    outMeshFilename = util::io::removeExtensionFromFile(util::io::getFilenameFromPath(meshFilename));
  }
  else {
    std::cerr << "testMPI() accepts either 1 or 2 for testOption. The value "
              << testOption << " is invalid.\n";
    exit(1);
  }
 
  // print mesh data and write mesh to a file
  util::io::print(mesh.printStr());
 
  // calculate nonlocal neighborhood
  double horizon = mHorizon*mesh.d_h;
  std::vector<std::vector<size_t>> nodeNeighs(mesh.d_numNodes);
  geom::computeNonlocalNeighborhood(mesh.d_nodes, horizon, nodeNeighs);
 
  // partition the mesh on root processor and broadcast to other processors
  mesh.d_nodePartition.resize(mesh.d_numNodes);
  util::io::print("\n\nCreating partition of mesh\n\n");
  if (mpiRank == 0)
    mesh::metisGraphPartition("metis_kway", &mesh, nodeNeighs, nPart);
 
  util::io::print("\n\nBroadcasting partition to all processors\n\n");
  MPI_Bcast(mesh.d_nodePartition.data(), mesh.d_numNodes,
            MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
 
 
  // Tasks:
  // 1. Create list of ghost nodes associated to neighboring processors
  // 2. Create a method that updates displacement of ghost nodes via MPI communication
  // 3. Create two set of nodes owned by this processor; internal set will have nodes
  // that do not depend on ghost nodes and boundary set that depend on ghost nodes
 
  // store id of nodes owned by this processor
  std::vector<size_t> ownedNodes, ownedInternalNodes, ownedBdryNodes;
 
  // for each neighboring processor, store id of ghost nodes owned by that processor
  std::vector<std::pair<std::vector<size_t>, std::vector<size_t>>> ghostData;
 
  // fill the owned and ghost node vectors
  util::io::print("\n\nCalling setupOwnerAndGhost()\n\n");
  setupOwnerAndGhost(mpiSize, mpiRank,
                     mesh.d_nodePartition, nodeNeighs,
                     ownedNodes, ownedInternalNodes, ownedBdryNodes,
                     ghostData);
 
  // create dummy displacement vector with random values
  std::vector<util::Point> dispNodes(mesh.d_numNodes, util::Point(-1., -1., -1.));
  {
    //  int seed = mpiRank;
    //  RandGenerator gen(util::get_rd_gen(seed));
    //  auto dist = util::DistributionSample<UniformDistribution>(0., 1., seed);
    //  for (auto i: ownedNodes)
    //    dispNodes[i] = util::Point(dist(), dist(), dist());
    // assign value to owned nodes
    for (auto i: ownedNodes)
      dispNodes[i] = util::Point(mpiRank + 1, (mpiRank + 1)*100, (mpiRank + 1)*10000);
  }
 
  // MPI communication to send and receive ghost nodes data
  printMsg("\n\nCalling exchangeDispData()\n\n", mpiRank, 0);
  std::vector<std::pair<std::vector<util::Point>, std::vector<util::Point>>> dispGhostData;
  exchangeDispData(mpiSize, mpiRank, ghostData, dispGhostData, dispNodes);
 
  if (true) {
    printMsg("\n\nDebugging dispGhostData()\n\n", mpiRank, -1);
    bool debug_failed = false;
    for (size_t j_proc = 0; j_proc < mpiSize; j_proc++) {
      auto &recvIds = ghostData[j_proc].first;
 
      if (j_proc != mpiRank and recvIds.size() != 0) {
        for (size_t k = 0; k < recvIds.size(); k++) {
          auto &uk = dispNodes[recvIds[k]];
 
          // verify we received correct value of uk
          if (uk[0] != j_proc + 1 or uk[1] != 100 * (j_proc + 1)
              or uk[2] != 10000 * (j_proc + 1)) {
            debug_failed = true;
            util::io::print(std::format("    MPI exchange error: j_proc = {}, "
                                     "uk = ({}, {}, {})\n", j_proc,
                                     uk[0], uk[1], uk[2]),
                     util::io::print_default_tab, -1);
          }
        }
      }
    } // loop over j_proc
 
    if (debug_failed)
      util::io::print(std::format("\n\nDEBUG failed for processor = {}\n\n", mpiRank), util::io::print_default_tab, -1);
    else
      util::io::print(std::format("\n\nDEBUG passed for processor = {}\n\n", mpiRank), util::io::print_default_tab, -1);
  }
}

References geom::computeNonlocalNeighborhood(), mesh::createUniformMesh(), anonymous_namespace{testParallelCompLib.cpp}::exchangeDispData(), util::io::getFilenameFromPath(), mesh::metisGraphPartition(), util::io::print(), util::io::print_default_tab, anonymous_namespace{testParallelCompLib.cpp}::printMsg(), util::io::removeExtensionFromFile(), and anonymous_namespace{testParallelCompLib.cpp}::setupOwnerAndGhost().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testNanoflann()

template<int dim = 3>

template std::string test::testNanoflann< 3 >	(	size_t	N,
		double	L,
		double	dL,
		int	seed
	)

Perform test on nsearch.

Parameters

N	size of particle cloud in each dimension. total size would be N^3
L	Size of unit cell to create crystal lattice point cloud
dL	Perturbation of lattice sites
seed	Seed
dim	Dimension

Returns

str Description of the search and its result

Definition at line 494 of file testNSearchLib.cpp.

                                                                     {
 
  if (dim < 2 or dim > 3) {
    return "testNanoflann: only dim = 2, 3 are accepted.\n";
  }
 
  // create 3D lattice and perturb each lattice point
  size_t Nx, Ny, Nz;
  Nx = Ny = Nz = N;
  size_t N_tot = Nx * Ny;
  if (dim == 3)
    N_tot = N_tot * Nz;
 
  std::vector<util::Point> x(N_tot, util::Point());
  lattice(L, Nx, Ny, Nz, dL, seed, x, dim);
  std::cout << "Total points = " << x.size() << "\n";
 
 
  if (dim == 2) {
    std::vector<std::vector<size_t>> neigh_nflann(N_tot, std::vector<size_t>());
    std::vector<std::vector<float>> neigh_nflann_sq_dist(N_tot,
                                                         std::vector<float>());
 
    std::vector<std::vector<size_t>> neigh_nflann_3d(N_tot, std::vector<size_t>());
    std::vector<std::vector<float>> neigh_nflann_3d_sq_dist(N_tot,
                                                         std::vector<float>());
 
    std::vector<std::vector<size_t>> neigh_brute(N_tot, std::vector<size_t>());
    std::vector<std::vector<float>> neigh_brute_sq_dist(N_tot,
                                                        std::vector<float>());
 
    // brute-force search
    double search_r = 1.5 * L;
    auto brute_force_search_time =
            neighSearchBrute(x, search_r, neigh_brute, neigh_brute_sq_dist);
 
    // nanoflann tree search
    std::unique_ptr<nsearch::NFlannSearchKd<3>> nflann_nsearch_3d
            = std::make_unique<nsearch::NFlannSearchKd<3>>(x, 0);
 
    std::unique_ptr<nsearch::NFlannSearchKd<dim>> nflann_nsearch
            = std::make_unique<nsearch::NFlannSearchKd<dim>>(x, 0);
 
    auto nflann_tree_set_time_3d = nflann_nsearch_3d->setInputCloud();
    auto nflann_tree_search_time_3d =
            neighSearchTreeSizet(x, nflann_nsearch_3d, search_r, neigh_nflann_3d, neigh_nflann_3d_sq_dist);
 
    auto nflann_tree_set_time = nflann_nsearch->setInputCloud();
    auto nflann_tree_search_time =
            neighSearchTreeSizet(x, nflann_nsearch, search_r, neigh_nflann, neigh_nflann_sq_dist);
 
    // Compare search results
    auto nflann_brute_compare = compare_results(
            neigh_nflann, neigh_brute, {"nflann_tree", "brute_force"}, -1, true);
 
    auto nflann_brute_compare_3d = compare_results(
            neigh_nflann_3d, neigh_brute, {"nflann_tree-3d", "brute_force"}, -1, true);
 
    std::ostringstream msg;
    msg << std::format("  Setup times (microseconds): \n"
                       "    nflann_tree_set_time = {} \n "
                       "    nflann_tree_set_time_3d = {}\n",
                       nflann_tree_set_time, nflann_tree_set_time_3d);
 
    msg << std::format("  Search times (microseconds): \n"
                       "    brute_force_search_time = {}\n"
                       "    nflann_tree_search_time = {}\n"
                       "    nflann_tree_search_time_3d = {}\n",
                       brute_force_search_time,
                       nflann_tree_search_time,
                       nflann_tree_search_time_3d);
 
    msg << std::format("  Comparison results: \n"
                       "    nflann_brute_compare: \n{}\n",
                       nflann_brute_compare);
 
    msg << std::format("  Comparison results: \n"
                       "    nflann_brute_compare_3d: \n{}\n",
                       nflann_brute_compare_3d);
 
    return msg.str();
  }
  else if (dim == 3) {
    std::vector<std::vector<size_t>> neigh_nflann(N_tot, std::vector<size_t>());
    std::vector<std::vector<float>> neigh_nflann_sq_dist(N_tot,
                                                         std::vector<float>());
 
    std::vector<std::vector<size_t>> neigh_brute(N_tot, std::vector<size_t>());
    std::vector<std::vector<float>> neigh_brute_sq_dist(N_tot,
                                                        std::vector<float>());
 
    // brute-force search
    double search_r = 1.5 * L;
    auto brute_force_search_time =
            neighSearchBrute(x, search_r, neigh_brute, neigh_brute_sq_dist);
 
    // nanoflann tree search
    std::unique_ptr<nsearch::NFlannSearchKd<dim>> nflann_nsearch
            = std::make_unique<nsearch::NFlannSearchKd<dim>>(x, 0);
 
    auto nflann_tree_set_time = nflann_nsearch->setInputCloud();
    auto nflann_tree_search_time =
            neighSearchTreeSizet(x, nflann_nsearch, search_r, neigh_nflann, neigh_nflann_sq_dist);
 
    // Compare search results
    auto nflann_brute_compare = compare_results(
            neigh_nflann, neigh_brute, {"nflann_tree", "brute_force"}, -1, true);
 
    std::ostringstream msg;
    msg << std::format("  Setup times (microseconds): \n"
                       "    nflann_tree_set_time = {}\n",
                       nflann_tree_set_time);
 
    msg << std::format("  Search times (microseconds): \n"
                       "    brute_force_search_time = {}\n"
                       "    nflann_tree_search_time = {}\n",
                       brute_force_search_time,
                       nflann_tree_search_time);
 
    msg << std::format("  Comparison results: \n"
                       "    nflann_brute_compare: \n{}\n",
                       nflann_brute_compare);
 
    return msg.str();
  }
 
 
}

References anonymous_namespace{testNSearchLib.cpp}::compare_results(), anonymous_namespace{testNSearchLib.cpp}::lattice(), anonymous_namespace{testNSearchLib.cpp}::neighSearchBrute(), and anonymous_namespace{testNSearchLib.cpp}::neighSearchTreeSizet().

Here is the call graph for this function:

testNanoflannClosestPoint()

std::string test::testNanoflannClosestPoint	(	size_t	N,
		double	L,
		double	dL,
		int	seed
	)

Perform test on nsearch.

Parameters

N	size of particle cloud in each dimension. total size would be N^3
L	Size of unit cell to create crystal lattice point cloud
dL	Perturbation of lattice sites
seed	Seed

Returns

str Description of the search and its result

Definition at line 781 of file testNSearchLib.cpp.

                                                                                 {
 
  // create 3D lattice and perturb each lattice point
  size_t Nx, Ny, Nz;
  Nx = Ny = Nz = N;
  size_t N_tot = Nx * Ny * Nz;
 
  std::vector<util::Point> x(N_tot, util::Point());
  lattice(L, Nx, Ny, Nz, dL, seed, x, 3);
  std::cout << "Total points = " << x.size() << "\n";
 
  // nanoflann tree search
  auto nflann_nsearch
          = std::make_unique<nsearch::NFlannSearchKd<3>>(x, 0);
 
  auto nflann_tree_set_time = nflann_nsearch->setInputCloud();
 
  std::vector<size_t> err_points;
  std::vector<util::Point> search_points;
  std::vector<double> err_dist;
  auto nsearch_search_time = neighSearchTreeClosestPointSizet(
          x, nflann_nsearch, seed, L, dL, search_points, err_points, err_dist);
 
    // Compare search results
    auto nflann_compare = compare_closest_point_results(x,
                                                              search_points,
                                                              err_points,
                                                              err_dist, false);
 
    std::ostringstream msg;
    msg << std::format("  Setup times (microseconds): \n"
                       "    nflann_tree_set_time = {}\n",
                       nflann_tree_set_time);
 
    msg << std::format("  Comparison results: \n"
                       "    nflann_compare: \n{}\n",
                       nflann_compare);
 
    return msg.str();
 
 
 
}

References anonymous_namespace{testNSearchLib.cpp}::compare_closest_point_results(), anonymous_namespace{testNSearchLib.cpp}::lattice(), and anonymous_namespace{testNSearchLib.cpp}::neighSearchTreeClosestPointSizet().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testNanoflannExcludeInclude()

template<int dim = 3>

template std::string test::testNanoflannExcludeInclude< 3 >	(	size_t	N,
		double	L,
		double	dL,
		int	seed,
		testNSearchData &	data
	)

Perform test on nsearch.

Parameters

N	size of particle cloud in each dimension. total size would be N^3
L	Size of unit cell to create crystal lattice point cloud
dL	Perturbation of lattice sites
seed	Seed
dim	Dimension
data	Search data

Returns

str Description of the search and its result

Definition at line 625 of file testNSearchLib.cpp.

                                                     {
 
  if (dim < 2 or dim > 3) {
    return "testNanoflannExcludeInclude: only dim = 2, 3 are accepted.\n";
  }
 
  // create 3D lattice and perturb each lattice point
  size_t Nx, Ny, Nz;
  Nx = Ny = Nz = N;
  size_t N_tot = Nx * Ny;
  if (dim == 3)
    N_tot = N_tot * Nz;
 
  std::vector<util::Point> x(N_tot, util::Point());
  std::vector<size_t> xTags(N_tot, 0);
  lattice(L, Nx, Ny, Nz, dL, seed, x, dim);
  assignRandomTags(x, data.d_numTags, seed, xTags);
  data.d_numPoints = x.size();
  std::cout << "Total points = " << x.size() << "\n";
 
  std::vector<std::vector<size_t>> neigh_default_nflann(N_tot, std::vector<size_t>());
  std::vector<std::vector<float>> neigh_default_nflann_sq_dist(N_tot,
                                                               std::vector<float>());
 
  std::vector<std::vector<size_t>> neigh_exclude_nflann(N_tot, std::vector<size_t>());
  std::vector<std::vector<float>> neigh_exclude_nflann_sq_dist(N_tot,
                                                       std::vector<float>());
 
  std::vector<std::vector<size_t>> neigh_include_nflann(N_tot, std::vector<size_t>());
  std::vector<std::vector<float>> neigh_include_nflann_sq_dist(N_tot,
                                                               std::vector<float>());
 
 
  std::vector<std::vector<size_t>> neigh_default_brute(N_tot, std::vector<size_t>());
  std::vector<std::vector<float>> neigh_default_brute_sq_dist(N_tot,
                                                              std::vector<float>());
 
  std::vector<std::vector<size_t>> neigh_exclude_brute(N_tot, std::vector<size_t>());
  std::vector<std::vector<float>> neigh_exclude_brute_sq_dist(N_tot,
                                                      std::vector<float>());
 
  std::vector<std::vector<size_t>> neigh_include_brute(N_tot, std::vector<size_t>());
  std::vector<std::vector<float>> neigh_include_brute_sq_dist(N_tot,
                                                              std::vector<float>());
 
  // brute-force search
  double search_r = 3. * L;
 
  data.d_defaultBruteSearchTime =
          neighSearchBrute(x,
                           search_r,
                           neigh_default_brute,
                           neigh_default_brute_sq_dist);
 
  data.d_excludeBruteSearchTime =
          neighSearchBruteExcludeInclude(x, xTags, search_r,
                                  neigh_exclude_brute,
                                  neigh_exclude_brute_sq_dist,
                                  1);
 
  data.d_includeBruteSearchTime =
          neighSearchBruteExcludeInclude(x, xTags, search_r,
                                  neigh_include_brute,
                                  neigh_include_brute_sq_dist,
                                  2);
 
  // nanoflann tree search
  std::unique_ptr<nsearch::NFlannSearchKd<dim>> nflann_nsearch
          = std::make_unique<nsearch::NFlannSearchKd<dim>>(x, 0,
                  data.d_leafMaxSize);
 
  data.d_treeBuildTime = nflann_nsearch->setInputCloud();
 
  // default
  data.d_defaultNFlannSearchTime =
          neighSearchTreeSizet(x, nflann_nsearch,
                               search_r, neigh_default_nflann,
                               neigh_default_nflann_sq_dist);
 
  // exclude
  data.d_excludeNFlannSearchTime =
          neighSearchTreeSizetExcludeInclude(x, xTags, nflann_nsearch,
                                      search_r, neigh_exclude_nflann,
                                      neigh_exclude_nflann_sq_dist,
                                      1);
 
  // include
  data.d_includeNFlannSearchTime =
          neighSearchTreeSizetExcludeInclude(x, xTags, nflann_nsearch,
                                      search_r, neigh_include_nflann,
                                      neigh_include_nflann_sq_dist,
                                      2);
 
  // Compare search results
  auto nflann_brute_compare_default = compare_results(
          neigh_default_nflann, neigh_default_brute,
          {"nflann_tree_default", "brute_force_default"}, -1, true);
 
  auto nflann_brute_compare_exclude = compare_results(
          neigh_exclude_nflann, neigh_exclude_brute,
          {"nflann_tree_exclude", "brute_force_exclude"}, -1, true);
 
  auto nflann_brute_compare_include = compare_results(
          neigh_include_nflann, neigh_include_brute,
          {"nflann_tree_include", "brute_force_include"}, -1, true);
 
  std::ostringstream msg;
  msg << std::format("  Setup times (microseconds): \n"
                     "    nflann_tree_set_time = {}\n",
                     data.d_treeBuildTime);
 
  msg << std::format("  Default search times (microseconds): \n"
                     "    brute_force_search_time = {}\n"
                     "    nflann_tree_search_time = {}\n",
                     data.d_defaultBruteSearchTime,
                     data.d_defaultNFlannSearchTime);
 
  msg << std::format("  Exclude comparison results: \n"
                     "    nflann_brute_compare: \n{}\n",
                     nflann_brute_compare_default);
 
  msg << std::format("  Exclude search times (microseconds): \n"
                     "    brute_force_search_time = {}\n"
                     "    nflann_tree_search_time = {}\n",
                     data.d_excludeBruteSearchTime,
                     data.d_excludeNFlannSearchTime);
 
  msg << std::format("  Exclude comparison results: \n"
                     "    nflann_brute_compare: \n{}\n",
                     nflann_brute_compare_exclude);
 
  msg << std::format("  Include search times (microseconds): \n"
                     "    brute_force_search_time = {}\n"
                     "    nflann_tree_search_time = {}\n",
                     data.d_includeBruteSearchTime,
                     data.d_includeNFlannSearchTime);
 
  msg << std::format("  Include comparison results: \n"
                     "    nflann_brute_compare: \n{}\n",
                     nflann_brute_compare_include);
 
 
  msg << std::format("  Nflann all search times (microseconds): \n"
                     "    default = {}\n"
                     "    exclude = {}\n"
                     "    include = {}\n",
                     data.d_defaultNFlannSearchTime,
                     data.d_excludeNFlannSearchTime,
                     data.d_includeNFlannSearchTime);
 
  return msg.str();
 
 
}

References anonymous_namespace{testNSearchLib.cpp}::assignRandomTags(), anonymous_namespace{testNSearchLib.cpp}::compare_results(), anonymous_namespace{testNSearchLib.cpp}::lattice(), anonymous_namespace{testNSearchLib.cpp}::neighSearchBrute(), anonymous_namespace{testNSearchLib.cpp}::neighSearchBruteExcludeInclude(), anonymous_namespace{testNSearchLib.cpp}::neighSearchTreeSizet(), and anonymous_namespace{testNSearchLib.cpp}::neighSearchTreeSizetExcludeInclude().

Here is the call graph for this function:

testPatchQuad()

void test::testPatchQuad

(

)

Perform test on quadrature points on line elements (NOT IMPLEMENTED)

This function performs accuracy test of the quadrature points for integration over reference line with vertices at {-1, 1}. List of tests are as follows:

1. Computes quadrature points of the given order, writes them to the file, and checks if the sum of quadrature weights is equal to 0.5 (area of reference triangle).

Also tests the exactness of the integration of the polynomial upto given order. Suppose \( n\) is the order of quadrature point, then we test if the integration of the function \( f(s,t) = s^\alpha\, t^\beta \) is exact for \( \alpha \) and \( \beta \) such that \( \alpha+\beta \leq n \). The exact integration of function \( f\) over reference triangle is

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds = \sum_{i=0}^{\beta+1} (-1)^i \frac{{{\beta + 1} \choose i}}{(\alpha + i +1) (\beta + 1)}, \]

where

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

We have \( {a \choose 0} = 1 \) so that term for \( i=0\) is not zero. Above formula gives the exact value of integral of \( f(s,t) = s^\alpha\, t^\beta \) over reference triangle. Approximation by quadrature point is as follows

\[ I_{approx} = \sum_{q=1}^{Q} w_q f(s_q, t_q) \]

where \(Q\) is the total number of quad points, \( w_q\) and \((s_q, t_q)\) are the \( q^{th} \) quad weight and point. In this test, we compare \( I_{exact} \) and \( I_{approx} \) and report problem if both do not match.

1. Test the accuracy for simple quadrangle mesh on square domain [0,1]^2. Exact integration of polynomial \( f(x,y) = x^\alpha \, y^\beta \) on \([0,1]^2\) is given by

   \[ I_{exact} = \frac{1}{(\alpha+1) (\beta+1)}. \]

   We compare above with the approximation computed from the quadrature points. We consider \(\alpha + \beta \leq n \) where \( n\) is the order of approximation we are testing.

Parameters

n	Order of quadrature point approximation
filepath	Path where mesh data for test can be found (expects files 'triMesh_nodes.csv' and 'triMesh_elements.csv' inside the filepath)

Definition at line 827 of file testFeLib.cpp.

                         {
  const double ptol = 1.0e-5;
  fe::QuadElem quad(1);
  const std::vector<util::Point> nodes = {
      util::Point(0., 0., 0.), util::Point(1., 0., 0.),
      util::Point(1., 1., 0.), util::Point(0., 1., 0.)};
  const std::vector<util::Point> u = {
      util::Point(0., 0., 0.), util::Point(2., 0., 0.),
      util::Point(2., 3., 0.), util::Point(0., 3., 0.)};
  auto qds = quad.getQuadDatas(nodes);
  size_t nfail = 0;
  for (const auto &qd : qds) {
    auto e = strainFromB(fe::B(qd.d_derShapes, 2), u);
    if (std::abs(e(0, 0) - 2.) > ptol || std::abs(e(1, 1) - 3.) > ptol ||
        std::abs(e(0, 1)) > ptol)
      nfail++;
  }
  std::cout << "**********************************\n";
  std::cout << "Patch test (quad, linear u)\n";
  std::cout << "**********************************\n";
  if (nfail) {
    std::cout << "PATCH QUAD : FAIL.\n";
    exit(EXIT_FAILURE);
  }
  std::cout << "PATCH QUAD : PASS.\n";
}

References fe::BaseElem::getQuadDatas(), and anonymous_namespace{testFeLib.cpp}::strainFromB().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testPatchTet()

void test::testPatchTet

(

)

Perform test on quadrature points on line elements (NOT IMPLEMENTED)

This function performs accuracy test of the quadrature points for integration over reference line with vertices at {-1, 1}. List of tests are as follows:

1. Computes quadrature points of the given order, writes them to the file, and checks if the sum of quadrature weights is equal to 0.5 (area of reference triangle).

Also tests the exactness of the integration of the polynomial upto given order. Suppose \( n\) is the order of quadrature point, then we test if the integration of the function \( f(s,t) = s^\alpha\, t^\beta \) is exact for \( \alpha \) and \( \beta \) such that \( \alpha+\beta \leq n \). The exact integration of function \( f\) over reference triangle is

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds = \sum_{i=0}^{\beta+1} (-1)^i \frac{{{\beta + 1} \choose i}}{(\alpha + i +1) (\beta + 1)}, \]

where

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

We have \( {a \choose 0} = 1 \) so that term for \( i=0\) is not zero. Above formula gives the exact value of integral of \( f(s,t) = s^\alpha\, t^\beta \) over reference triangle. Approximation by quadrature point is as follows

\[ I_{approx} = \sum_{q=1}^{Q} w_q f(s_q, t_q) \]

where \(Q\) is the total number of quad points, \( w_q\) and \((s_q, t_q)\) are the \( q^{th} \) quad weight and point. In this test, we compare \( I_{exact} \) and \( I_{approx} \) and report problem if both do not match.

1. Test the accuracy for simple quadrangle mesh on square domain [0,1]^2. Exact integration of polynomial \( f(x,y) = x^\alpha \, y^\beta \) on \([0,1]^2\) is given by

   \[ I_{exact} = \frac{1}{(\alpha+1) (\beta+1)}. \]

   We compare above with the approximation computed from the quadrature points. We consider \(\alpha + \beta \leq n \) where \( n\) is the order of approximation we are testing.

Parameters

n	Order of quadrature point approximation
filepath	Path where mesh data for test can be found (expects files 'triMesh_nodes.csv' and 'triMesh_elements.csv' inside the filepath)

Definition at line 854 of file testFeLib.cpp.

                        {
  const double ptol = 1.0e-5;
  fe::TetElem tet(1);
  const std::vector<util::Point> nodes = {
      util::Point(0., 0., 0.), util::Point(1., 0., 0.),
      util::Point(0., 1., 0.), util::Point(0., 0., 1.)};
  const std::vector<util::Point> u = {
      util::Point(0., 0., 0.), util::Point(2., 0., 0.),
      util::Point(0., 3., 0.), util::Point(0., 0., 4.)};
  auto qds = tet.getQuadDatas(nodes);
  size_t nfail = 0;
  for (const auto &qd : qds) {
    auto e = strainFromB(fe::B(qd.d_derShapes, 3), u);
    if (std::abs(e(0, 0) - 2.) > ptol || std::abs(e(1, 1) - 3.) > ptol ||
        std::abs(e(2, 2) - 4.) > ptol || std::abs(e(0, 1)) > ptol ||
        std::abs(e(0, 2)) > ptol || std::abs(e(1, 2)) > ptol)
      nfail++;
  }
  std::cout << "**********************************\n";
  std::cout << "Patch test (tet, linear u)\n";
  std::cout << "**********************************\n";
  if (nfail) {
    std::cout << "PATCH TET : FAIL.\n";
    exit(EXIT_FAILURE);
  }
  std::cout << "PATCH TET : PASS.\n";
}

References fe::BaseElem::getQuadDatas(), and anonymous_namespace{testFeLib.cpp}::strainFromB().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testPatchTri()

void test::testPatchTri

(

)

Perform test on quadrature points on line elements (NOT IMPLEMENTED)

This function performs accuracy test of the quadrature points for integration over reference line with vertices at {-1, 1}. List of tests are as follows:

1. Computes quadrature points of the given order, writes them to the file, and checks if the sum of quadrature weights is equal to 0.5 (area of reference triangle).

Also tests the exactness of the integration of the polynomial upto given order. Suppose \( n\) is the order of quadrature point, then we test if the integration of the function \( f(s,t) = s^\alpha\, t^\beta \) is exact for \( \alpha \) and \( \beta \) such that \( \alpha+\beta \leq n \). The exact integration of function \( f\) over reference triangle is

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds = \sum_{i=0}^{\beta+1} (-1)^i \frac{{{\beta + 1} \choose i}}{(\alpha + i +1) (\beta + 1)}, \]

where

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

We have \( {a \choose 0} = 1 \) so that term for \( i=0\) is not zero. Above formula gives the exact value of integral of \( f(s,t) = s^\alpha\, t^\beta \) over reference triangle. Approximation by quadrature point is as follows

\[ I_{approx} = \sum_{q=1}^{Q} w_q f(s_q, t_q) \]

where \(Q\) is the total number of quad points, \( w_q\) and \((s_q, t_q)\) are the \( q^{th} \) quad weight and point. In this test, we compare \( I_{exact} \) and \( I_{approx} \) and report problem if both do not match.

1. Test the accuracy for simple quadrangle mesh on square domain [0,1]^2. Exact integration of polynomial \( f(x,y) = x^\alpha \, y^\beta \) on \([0,1]^2\) is given by

   \[ I_{exact} = \frac{1}{(\alpha+1) (\beta+1)}. \]

   We compare above with the approximation computed from the quadrature points. We consider \(\alpha + \beta \leq n \) where \( n\) is the order of approximation we are testing.

Parameters

n	Order of quadrature point approximation
filepath	Path where mesh data for test can be found (expects files 'triMesh_nodes.csv' and 'triMesh_elements.csv' inside the filepath)

Definition at line 800 of file testFeLib.cpp.

                        {
  const double ptol = 1.0e-5;
  fe::TriElem tri(1);
  const std::vector<util::Point> nodes = {
      util::Point(0., 0., 0.), util::Point(1., 0., 0.),
      util::Point(0., 1., 0.)};
  const std::vector<util::Point> u = {
      util::Point(0., 0., 0.), util::Point(2., 0., 0.),
      util::Point(0., 3., 0.)};
  auto qds = tri.getQuadDatas(nodes);
  size_t nfail = 0;
  for (const auto &qd : qds) {
    auto e = strainFromB(fe::B(qd.d_derShapes, 2), u);
    if (std::abs(e(0, 0) - 2.) > ptol || std::abs(e(1, 1) - 3.) > ptol ||
        std::abs(e(0, 1)) > ptol)
      nfail++;
  }
  std::cout << "**********************************\n";
  std::cout << "Patch test (triangle, linear u)\n";
  std::cout << "**********************************\n";
  if (nfail) {
    std::cout << "PATCH TRI : FAIL.\n";
    exit(EXIT_FAILURE);
  }
  std::cout << "PATCH TRI : PASS.\n";
}

References fe::BaseElem::getQuadDatas(), and anonymous_namespace{testFeLib.cpp}::strainFromB().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testPatchTriDistorted()

void test::testPatchTriDistorted

(

)

Perform test on quadrature points on line elements (NOT IMPLEMENTED)

This function performs accuracy test of the quadrature points for integration over reference line with vertices at {-1, 1}. List of tests are as follows:

1. Computes quadrature points of the given order, writes them to the file, and checks if the sum of quadrature weights is equal to 0.5 (area of reference triangle).

Also tests the exactness of the integration of the polynomial upto given order. Suppose \( n\) is the order of quadrature point, then we test if the integration of the function \( f(s,t) = s^\alpha\, t^\beta \) is exact for \( \alpha \) and \( \beta \) such that \( \alpha+\beta \leq n \). The exact integration of function \( f\) over reference triangle is

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds = \sum_{i=0}^{\beta+1} (-1)^i \frac{{{\beta + 1} \choose i}}{(\alpha + i +1) (\beta + 1)}, \]

where

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

We have \( {a \choose 0} = 1 \) so that term for \( i=0\) is not zero. Above formula gives the exact value of integral of \( f(s,t) = s^\alpha\, t^\beta \) over reference triangle. Approximation by quadrature point is as follows

\[ I_{approx} = \sum_{q=1}^{Q} w_q f(s_q, t_q) \]

where \(Q\) is the total number of quad points, \( w_q\) and \((s_q, t_q)\) are the \( q^{th} \) quad weight and point. In this test, we compare \( I_{exact} \) and \( I_{approx} \) and report problem if both do not match.

1. Test the accuracy for simple quadrangle mesh on square domain [0,1]^2. Exact integration of polynomial \( f(x,y) = x^\alpha \, y^\beta \) on \([0,1]^2\) is given by

   \[ I_{exact} = \frac{1}{(\alpha+1) (\beta+1)}. \]

   We compare above with the approximation computed from the quadrature points. We consider \(\alpha + \beta \leq n \) where \( n\) is the order of approximation we are testing.

Parameters

n	Order of quadrature point approximation
filepath	Path where mesh data for test can be found (expects files 'triMesh_nodes.csv' and 'triMesh_elements.csv' inside the filepath)

Definition at line 882 of file testFeLib.cpp.

                                 {
  const double ptol = 1.0e-5;
  fe::TriElem tri(1);
  const std::vector<util::Point> nodes = {
      util::Point(0., 0., 0.), util::Point(1.3, 0.2, 0.),
      util::Point(0.15, 1.1, 0.)};
  const std::vector<util::Point> u = {
      util::Point(0., 0., 0.), util::Point(2.6, 0.6, 0.),
      util::Point(0.3, 3.3, 0.)};
  auto qds = tri.getQuadDatas(nodes);
  size_t nfail = 0;
  for (const auto &qd : qds) {
    auto e = strainFromB(fe::B(qd.d_derShapes, 2), u);
    if (std::abs(e(0, 0) - 2.) > ptol || std::abs(e(1, 1) - 3.) > ptol ||
        std::abs(e(0, 1)) > ptol)
      nfail++;
  }
  std::cout << "**********************************\n";
  std::cout << "Patch test (distorted triangle, linear u)\n";
  std::cout << "**********************************\n";
  if (nfail) {
    std::cout << "PATCH TRI DISTORTED : FAIL.\n";
    exit(EXIT_FAILURE);
  }
  std::cout << "PATCH TRI DISTORTED : PASS.\n";
}

References fe::BaseElem::getQuadDatas(), and anonymous_namespace{testFeLib.cpp}::strainFromB().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testQuadElem()

void test::testQuadElem	(	size_t	n,
		std::string	filepath
	)

Perform test on quadrature points on quadrangle elements.

This function performs accuracy test of the quadrature points for integration over reference triangle with vertices at {(0,0), (1,0), (0,1)}. List of tests are as follows:

1. Computes quadrature points of the given order, writes them to the file, and checks if the sum of quadrature weights is equal to 2, i.e. area of reference quadrangle element.

Also tests the exactness of the integration of the polynomial upto given order. Suppose \( n\) is the order of quadrature point, then we test if the integration of the function \( f(s,t) = s^\alpha\, t^\beta \) is exact for \( \alpha \) and \( \beta \) such that \( \alpha+\beta \leq n \). The exact integration of function \( f\) over reference triangle is

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds = \sum_{i=0}^{\beta+1} (-1)^i \frac{{{\beta + 1} \choose i}}{(\alpha + i +1) (\beta + 1)}, \]

where

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

We have \( {a \choose 0} = 1 \) so that term for \( i=0\) is not zero. Above formula gives the exact value of integral of \( f(s,t) = s^\alpha\, t^\beta \) over reference triangle. Approximation by quadrature point is as follows

\[ I_{approx} = \sum_{q=1}^{Q} w_q f(s_q, t_q) \]

where \(Q\) is the total number of quad points, \( w_q\) and \((s_q, t_q)\) are the \( q^{th} \) quad weight and point. In this test, we compare \( I_{exact} \) and \( I_{approx} \) and report problem if both do not match.

1. Test the accuracy for simple quadrangle mesh on square domain [0,1]^2. Exact integration of polynomial \( f(x,y) = x^\alpha \, y^\beta \) on \([0,1]^2\) is given by

   \[ I_{exact} = \frac{1}{(\alpha+1) (\beta+1)}. \]

   We compare above with the approximation computed from the quadrature points. We consider \(\alpha + \beta \leq n \) where \( n\) is the order of approximation we are testing.

Parameters

n	Order of quadrature point approximation
filepath	Path where mesh data for test can be found (expects files 'quadMesh_nodes.csv' and 'quadMesh_elements.csv' inside the filepath)

Definition at line 384 of file testFeLib.cpp.

                                                    {
 
  //
  // Test1: We test accuracy of integrals of polynomials over reference
  // quadrangle. Reference triangle {(-1,-1), (1,-1), (1,1), (-1,1)}.
  //
  // Test2: We consider simple mesh in meshFeTest.txt over square domain
  // [0,1]^2 and test the accuracy of polynomials over square domain.
  //
 
  // get Quadrature
  auto quad = fe::QuadElem(n);
 
  //
  // Test 1
  //
  size_t error_test_1 = 0;
  {
    // T1 (reference quadrangle)
    // get quad points at reference triangle
    std::vector<util::Point> nodes = {
        util::Point(-1., -1., 0.), util::Point(1., -1., 0.),
        util::Point(1., 1., 0.), util::Point(-1., 1., 0.)};
    std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
    double sum = 0.;
    for (auto qd : qds)
      sum += qd.d_w;
 
    if (std::abs(sum - 4.0) > tol) {
      std::cout << "Error in order = " << n
                << ". Sum of quad weights is not "
                   "equal to area of reference "
                   "quadrangle.\n";
      error_test_1++;
    }
 
    //
    // test the exactness of integration for polynomial
    //
    for (size_t i = 0; i <= 2 * n - 1; i++)
      for (size_t j = 0; j <= 2 * n - 1; j++) {
 
        //
        // when {(-1,-1), (1,-1), (1,1), (-1,1)}
        //
        nodes = {util::Point(-1., -1., 0.), util::Point(1., -1., 0.),
                 util::Point(1., 1., 0.), util::Point(-1., 1., 0.)};
        qds = quad.getQuadPoints(nodes);
        // test integration of polynomial f(s,t) = s^i t^j
        // get the exact integration
        double I_exact = test::getExactIntegrationRefQuad(i, j);
        if (!checkRefIntegration(n, i, j, qds, I_exact))
          error_test_1++;
 
        //
        // when {(-1,1), (-1,-1), (1,-1), (1,1)}
        //
        nodes = {util::Point(-1., 1., 0.), util::Point(-1., -1., 0.),
                 util::Point(1., -1., 0.), util::Point(1., 1., 0.)};
        qds = quad.getQuadPoints(nodes);
        //
        // After changing the order of vertices, we have got a new
        // triangle which is in coordinate system (x,y) and we are
        // integrating function f(x,y) = x^i y^j
        //
        // The quad data we have got is such that quad point is in (x,y)
        // coordinate, weight is such that determinant of the Jacobian is
        // included in the weight.
        //
        // Thus the following method for I_approx is correct.
        if (!checkRefIntegration(n, i, j, qds, I_exact))
          error_test_1++;
 
        //
        // when {(1,1), (-1,1), (-1,-1), (1,-1)}
        //
        nodes = {util::Point(1., 1., 0.), util::Point(-1., 1., 0.),
                 util::Point(-1., -1., 0.), util::Point(1., -1., 0.)};
        qds = quad.getQuadPoints(nodes);
        if (!checkRefIntegration(n, i, j, qds, I_exact))
          error_test_1++;
 
        //
        // when {(1,-1), (1,1), (-1,1), (-1,-1)}
        //
        nodes = {util::Point(1., -1., 0.), util::Point(1., 1., 0.),
                 util::Point(-1., 1., 0.), util::Point(-1., -1., 0.)};
        qds = quad.getQuadPoints(nodes);
        if (!checkRefIntegration(n, i, j, qds, I_exact))
          error_test_1++;
      }
  } // Test 1
 
  //
  // Test 2
  //
  size_t error_test_2 = 0;
  {
    static std::vector<util::Point> nodes;
    static std::vector<size_t> elements;
    static size_t num_vertex = 4;
    static size_t elem_type = util::vtk_type_quad;
    static size_t num_elems = 0;
    if (num_elems == 0) {
      readNodes(filepath + "/quadMesh_nodes.csv", nodes);
      num_elems = readElements(filepath + "/quadMesh_elements.csv", elem_type,
                               elements);
    }
 
    // loop over polynomials
    for (size_t i = 0; i <= 2 * n - 1; i++)
      for (size_t j = 0; j <= 2 * n - 1; j++) {
 
        double I_exact = 1. / (double(i + 1) * double(j + 1));
        double I_approx = 0.;
        // loop over elements and compute I_approx
        for (size_t e = 0; e < num_elems; e++) {
          std::vector<util::Point> enodes = {
              nodes[elements[num_vertex * e + 0]],
              nodes[elements[num_vertex * e + 1]],
              nodes[elements[num_vertex * e + 2]],
              nodes[elements[num_vertex * e + 3]]};
          std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
          for (auto qd : qds)
            I_approx +=
                qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j);
        }
 
        if (std::abs(I_exact - I_approx) > tol) {
          std::cout << "Error in order = " << n
                    << ". Exact integration = " << I_exact
                    << " and approximate integration = " << I_approx
                    << " of polynomial of order (i = " << i << " + j = " << j
                    << ") = " << i + j << " over square domain [0,1]x[0,1] "
                    << "is not matching using quadrature points.\n";
 
          error_test_2++;
        }
      }
  }
 
  if (n == 1) {
    std::cout << "**********************************\n";
    std::cout << "Quadrangle Quadrature Test\n";
    std::cout << "**********************************\n";
  }
  std::cout << "Quad order = " << n << ". ";
  std::cout << (error_test_1 == 0 ? "TEST 1 : PASS. " : "TEST 1 : FAIL. ");
  std::cout << (error_test_2 == 0 ? "TEST 2 : PASS. \n" : "TEST 2 : FAIL. \n");
}

References anonymous_namespace{testFeLib.cpp}::checkRefIntegration(), getExactIntegrationRefQuad(), anonymous_namespace{testFeLib.cpp}::readElements(), anonymous_namespace{testFeLib.cpp}::readNodes(), anonymous_namespace{testFeLib.cpp}::tol, and util::vtk_type_quad.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testTaskflow()

std::string test::testTaskflow	(	size_t	N,
		int	seed
	)

Perform test on taskflow.

Parameters

N	size of vector to profile taskflow
seed	Seed

Returns

str Description of the test and its result

Definition at line 221 of file testParallelCompLib.cpp.

                                               {
 
  auto nThreads = util::parallel::getNThreads();
  util::io::print(std::format("\n\ntestTaskflow(): Number of threads = {}\n\n", nThreads));
 
  // task: perform N computations in serial and using taskflow for_each
 
  // generate vector of random numbers
  RandGenerator gen(util::get_rd_gen(seed));
  auto dist = util::DistributionSample<UniformDistribution>(0., 1., seed);
 
  std::vector<double> x(N);
  std::vector<double> y1(N);
  std::vector<double> y2(N);
  std::generate(std::begin(x), std::end(x), [&] { return dist(); });
 
  // now do serial calculation
  auto t1 = steady_clock::now();
  for (size_t i = 0; i < N; i++) {
    if (x[i] < 0.5)
      y1[i] = f1(x[i]);
    else
      y1[i] = f2(x[i]);
  }
  auto t2 = steady_clock::now();
  auto dt12 = util::methods::timeDiff(t1, t2, "microseconds");
 
  // now do parallel calculation using taskflow
  tf::Executor executor(nThreads);
  tf::Taskflow taskflow;
 
  taskflow.for_each_index((std::size_t) 0, N, (std::size_t) 1, [&x, &y2](std::size_t i) {
      if (x[i] < 0.5)
        y2[i] = f1(x[i]);
      else
        y2[i] = f2(x[i]);
                          }
  ); // for_each
 
  executor.run(taskflow).get();
  auto t3 = steady_clock::now();
  auto dt23 = util::methods::timeDiff(t2, t3, "microseconds");
 
  // compare results
  double y_err = 0.;
  for (size_t i=0; i<N; i++)
    y_err += std::pow(y1[i] - y2[i], 2);
 
  if (y_err > 1.e-10) {
    std::cerr << std::format("Error: Serial and taskflow computation results do not match (squared error = {})\n",
                            y_err);
    exit(1);
  }
 
  // get time
  std::ostringstream msg;
  msg << std::format("  Serial computation took = {}ms\n", dt12);
  msg << std::format("  Taskflow computation took = {}ms\n", dt23);
  msg << std::format("  Speed-up factor = {}\n\n\n", dt12/dt23);
 
  return msg.str();
}

References anonymous_namespace{testParallelCompLib.cpp}::f1(), anonymous_namespace{testParallelCompLib.cpp}::f2(), util::get_rd_gen(), util::parallel::getNThreads(), util::io::print(), and util::methods::timeDiff().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testTetElem()

void test::testTetElem	(	size_t	n,
		std::string	filepath
	)

Perform test on quadrature points on tetrahedral elements.

Parameters

n	Order of quadrature point approximation
filepath	Path where mesh data for test can be found (expects files 'tetMesh_nodes.csv' and 'tetMesh_elements.csv' inside the filepath)

Definition at line 608 of file testFeLib.cpp.

                                                   {
 
  //
  // Test1: We test accuracy of integrals of polynomials over reference
  // tetrahedron. Reference element {(0,0,0), (1,0,0), (0,1,0), (0,0,1)}.
  //
  // Test2: We consider simple mesh in meshFeTest.txt over cubic domain
  // [0,1]^3 and test the accuracy of polynomials over cubic domain.
  //
 
  // get Quadrature
  auto quad = fe::TetElem(n);
 
  //
  // Test 1
  //
  size_t error_test_1 = 0;
  {
    // T1 (reference triangle)
    // get quad points at reference triangle
    std::vector<util::Point> nodes = {util::Point(), util::Point(1., 0., 0.),
                                       util::Point(0., 1., 0.),
                                       util::Point(0., 0., 1.)};
    std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
 
    double sum = 0.;
    for (auto qd : qds)
      sum += qd.d_w;
 
    if (std::abs(sum - 1. / 6.) > tol) {
      std::cout << "Error in order = " << n
                << ". Sum of quad weights is not "
                   "equal to volume of reference "
                   "tetrahedron.\n";
      error_test_1++;
    }
 
    //
    // test the exactness of integration for polynomial
    //
    for (size_t i = 0; i <= n; i++)
      for (size_t j = 0; j <= n; j++)
        for (size_t k = 0; k <= n; k++) {
 
          if (i + j + k > n)
            continue;
 
          //
          // +ve order of indices are:
          // {0,1,2,3}; {1,2,0,3}; {2,3,0,1}; {0,3,1,2}
 
          //
          // when {(0,0,0), (1,0,0), (0,1,0), (0,0,1)}
          //
          nodes = {util::Point(), util::Point(1., 0., 0.),
                   util::Point(0., 1., 0.), util::Point(0., 0., 1.)};
          qds = quad.getQuadPoints(nodes);
          // test integration of polynomial f(s,t) = s^i t^j
          // get the exact integration
          debug_id = 0;
          double I_exact = test::getExactIntegrationRefTet(i, j, k);
          if (!checkRefIntegration(n, i, j, k, qds, I_exact)) {
            error_test_1++;
          }
 
          //
          // when vertices are {(1,0,0), (0,1,0), (0,0,1), (0,0,0)}
          //
          nodes = {util::Point(1., 0., 0.), util::Point(0., 1., 0.),
                   util::Point(0., 0., 0.), util::Point(0., 0., 1.)};
          qds = quad.getQuadPoints(nodes);
          //
          // After changing the order of vertices, we have got a new
          // triangle which is in coordinate system (x,y) and we are
          // integrating function f(x,y) = x^i y^j
          //
          // The quad data we have got is such that quad point is in (x,y)
          // coordinate, weight is such that determinant of the Jacobian is
          // included in the weight.
          //
          // Thus the following method for I_approx is correct.
          debug_id = 1;
          if (!checkRefIntegration(n, i, j, k, qds, I_exact))
            error_test_1++;
 
          //
          // when vertices are {(0,1,0), (0,0,1), (0,0,0), (1,0,0)}
          //
          nodes = {util::Point(0., 1., 0.), util::Point(0., 0., 1.),
                   util::Point(0., 0., 0.), util::Point(1., 0., 0.)};
          qds = quad.getQuadPoints(nodes);
          debug_id = 2;
          if (!checkRefIntegration(n, i, j, k, qds, I_exact))
            error_test_1++;
 
          //
          // when vertices are {(0,0,1), (0,0,0), (1,0,0), (0,1,0)}
          //
          nodes = {util::Point(0., 0., 0.), util::Point(0., 0., 1.),
                   util::Point(1., 0., 0.), util::Point(0., 1., 0.)};
          qds = quad.getQuadPoints(nodes);
          debug_id = 3;
          if (!checkRefIntegration(n, i, j, k, qds, I_exact))
            error_test_1++;
        }
  } // Test 1
 
  //
  // Test 2
  //
  size_t error_test_2 = 0;
  if (false) {
    static std::vector<util::Point> nodes;
    static std::vector<size_t> elements;
    static size_t num_vertex = 4;
    static size_t elem_type = util::vtk_type_tetra;
    static size_t num_elems = 0;
    if (num_elems == 0) {
      readNodes(filepath + "tetMesh_nodes.csv", nodes);
      num_elems = readElements(filepath + "tetMesh_elements.csv", elem_type,
                               elements);
    }
 
    // loop over polynomials
    for (size_t i = 0; i <= n; i++)
      for (size_t j = 0; j <= n; j++)
        for (size_t k = 0; k <= n; k++) {
 
          if (i + j + k > n)
            continue;
 
          double I_exact = 1. / (double(i + 1) * double(j + 1) * double(k + 1));
          double I_approx = 0.;
          // loop over elements and compute I_approx
          for (size_t e = 0; e < num_elems; e++) {
            std::vector<util::Point> enodes = {
                nodes[elements[num_vertex * e + 0]],
                nodes[elements[num_vertex * e + 1]],
                nodes[elements[num_vertex * e + 2]],
                nodes[elements[num_vertex * e + 3]]};
            std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
            for (auto qd : qds) {
              I_approx += qd.d_w * std::pow(qd.d_p.d_x, i) *
                          std::pow(qd.d_p.d_y, j) * std::pow(qd.d_p.d_z, k);
 
              if (false) {
 
                std::cout << "Print " << i << " " << j << " " << k << "\n";
                std::cout << util::io::printStr(enodes) << "\n";
                std::vector<size_t> enode_ids = {
                      elements[num_vertex * e + 0],
                      elements[num_vertex * e + 1],
                      elements[num_vertex * e + 2],
                      elements[num_vertex * e + 3]};
                std::cout << util::io::printStr(enode_ids) << "\n";
                std::cout << qd.printStr() << "\n";
              }
 
            }
          }
 
          if (std::abs(I_exact - I_approx) > tol) {
            std::cout << "Error in order = " << n
                      << ". Exact integration = " << I_exact
                      << " and approximate integration = " << I_approx
                      << " of polynomial of order (i = " << i << " + j = " << j
                      << " + k = " << k << ") = " << i + j + k
                      << " over cubic domain [0,1]x[0,1]x[0,1] "
                      << "is not matching using quadrature points.\n";
 
            error_test_2++;
          }
        }
  }
 
  if (n == 1) {
    std::cout << "**********************************\n";
    std::cout << "Tetrahedron Quadrature Test\n";
    std::cout << "**********************************\n";
  }
  std::cout << "Quad order = " << n << ". ";
  if (error_test_1 == 0)
    std::cout << "TEST 1 : PASS. ";
  else
    std::cout << "TEST 1 : FAIL. ";
  //  if (error_test_2 == 0)
  //    std::cout << "TEST 2 : PASS. ";
  //  else
  //    std::cout << "TEST 2 : FAIL. ";
  std::cout << "\n";
}

References anonymous_namespace{testFeLib.cpp}::checkRefIntegration(), anonymous_namespace{testFeLib.cpp}::debug_id, getExactIntegrationRefTet(), util::io::printStr(), anonymous_namespace{testFeLib.cpp}::readElements(), anonymous_namespace{testFeLib.cpp}::readNodes(), anonymous_namespace{testFeLib.cpp}::tol, and util::vtk_type_tetra.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testTriElem()

void test::testTriElem	(	size_t	n,
		std::string	filepath
	)

Perform test on quadrature points on triangle elements.

This function performs accuracy test of the quadrature points for integration over reference triangle with vertices at {(0,0), (1,0), (0,1)}. List of tests are as follows:

1. Computes quadrature points of the given order, writes them to the file, and checks if the sum of quadrature weights is equal to 0.5, i.e. area of reference triangle.

Also tests the exactness of the integration of the polynomial upto given order. Suppose \( n\) is the order of quadrature point, then we test if the integration of the function \( f(s,t) = s^\alpha\, t^\beta \) is exact for \( \alpha \) and \( \beta \) such that \( \alpha+\beta \leq n \). The exact integration of function \( f\) over reference triangle is

\[ I_{exact} = \int_0^1 \int_0^{1-s} s^\alpha\, t^\beta \, dt\, ds = \sum_{i=0}^{\beta+1} (-1)^i \frac{{{\beta + 1} \choose i}}{(\alpha + i +1) (\beta + 1)}, \]

where

\[ {a \choose b} = \frac{a (a-1) (a-2) ... (a-b+1)}{1*2*3 ... *b}. \]

We have \( {a \choose 0} = 1 \) so that term for \( i=0\) is not zero. Above formula gives the exact value of integral of \( f(s,t) = s^\alpha\, t^\beta \) over reference triangle. Approximation by quadrature point is as follows

\[ I_{approx} = \sum_{q=1}^{Q} w_q f(s_q, t_q) \]

where \(Q\) is the total number of quad points, \( w_q\) and \((s_q, t_q)\) are the \( q^{th} \) quad weight and point. In this test, we compare \( I_{exact} \) and \( I_{approx} \) and report problem if both do not match.

1. Test the accuracy for simple triangular mesh on square domain [0,1]^2. Exact integration of polynomial \( f(x,y) = x^\alpha \, y^\beta \) on \([0,1]^2\) is given by

   \[ I_{exact} = \frac{1}{(\alpha+1) (\beta+1)}. \]

   We compare above with the approximation computed from the quadrature points. We consider \(\alpha + \beta \leq n \) where \( n\) is the order of approximation we are testing.

Parameters

n	Order of quadrature point approximation
filepath	Path where mesh data for test can be found (expects files 'triMesh_nodes.csv' and 'triMesh_elements.csv' inside the filepath)

Definition at line 231 of file testFeLib.cpp.

                                                   {
 
  //
  // Test1: We test accuracy of integrals of polynomials over reference
  // triangle. Reference triangle {(0,0), (1,0), (0,1)}.
  //
  // Test2: We consider simple mesh in meshFeTest.txt over square domain
  // [0,1]^2 and test the accuracy of polynomials over square domain.
  //
 
  // get Quadrature
  auto quad = fe::TriElem(n);
 
  //
  // Test 1
  //
  size_t error_test_1 = 0;
  {
    // T1 (reference triangle)
    // get quad points at reference triangle
    std::vector<util::Point> nodes = {util::Point(), util::Point(1., 0., 0.),
                                       util::Point(0., 1., 0.)};
    std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
    double sum = 0.;
    for (auto qd : qds)
      sum += qd.d_w;
 
    if (std::abs(sum - 0.5) > tol) {
      std::cout << "Error in order = " << n
                << ". Sum of quad weights is not "
                   "equal to area of reference "
                   "triangle.\n";
      error_test_1++;
    }
 
    //
    // test the exactness of integration for polynomial
    //
    for (size_t i = 0; i <= n; i++)
      for (size_t j = 0; j <= n; j++) {
 
        if (i + j > n)
          continue;
 
        //
        // when {(0,0), (1,0), (0,1)}
        //
        nodes = {util::Point(), util::Point(1., 0., 0.),
                 util::Point(0., 1., 0.)};
        qds = quad.getQuadPoints(nodes);
        // test integration of polynomial f(s,t) = s^i t^j
        // get the exact integration
        double I_exact = test::getExactIntegrationRefTri(i, j);
        if (!checkRefIntegration(n, i, j, qds, I_exact))
          error_test_1++;
 
        //
        // when vertices are {(1,0), (0,1), (0,0)}
        //
        nodes = {util::Point(1., 0., 0.), util::Point(0., 1., 0.),
                 util::Point()};
        qds = quad.getQuadPoints(nodes);
        //
        // After changing the order of vertices, we have got a new
        // triangle which is in coordinate system (x,y) and we are
        // integrating function f(x,y) = x^i y^j
        //
        // The quad data we have got is such that quad point is in (x,y)
        // coordinate, weight is such that determinant of the Jacobian is
        // included in the weight.
        //
        // Thus the following method for I_approx is correct.
        if (!checkRefIntegration(n, i, j, qds, I_exact))
          error_test_1++;
 
        //
        // when vertices are {(0,1), (0,0), (1,0)}
        //
        nodes = {util::Point(0., 1., 0.), util::Point(),
                 util::Point(1., 0., 0.)};
        qds = quad.getQuadPoints(nodes);
        if (!checkRefIntegration(n, i, j, qds, I_exact))
          error_test_1++;
      }
  } // Test 1
 
  //
  // Test 2
  //
  size_t error_test_2 = 0;
  {
    static std::vector<util::Point> nodes;
    static std::vector<size_t> elements;
    static size_t num_vertex = 3;
    static size_t elem_type = util::vtk_type_triangle;
    static size_t num_elems = 0;
    if (num_elems == 0) {
      readNodes(filepath + "/triMesh_nodes.csv", nodes);
      num_elems = readElements(filepath + "/triMesh_elements.csv", elem_type,
                               elements);
    }
 
    // loop over polynomials
    for (size_t i = 0; i <= n; i++)
      for (size_t j = 0; j <= n; j++) {
 
        if (i + j > n)
          continue;
 
        double I_exact = 1. / (double(i + 1) * double(j + 1));
        double I_approx = 0.;
        // loop over elements and compute I_approx
        for (size_t e = 0; e < num_elems; e++) {
          std::vector<util::Point> enodes = {
              nodes[elements[num_vertex * e + 0]],
              nodes[elements[num_vertex * e + 1]],
              nodes[elements[num_vertex * e + 2]]};
          std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
          for (auto qd : qds)
            I_approx +=
                qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j);
        }
 
        if (std::abs(I_exact - I_approx) > tol) {
          std::cout << "Error in order = " << n
                    << ". Exact integration = " << I_exact
                    << " and approximate integration = " << I_approx
                    << " of polynomial of order (i = " << i << " + j = " << j
                    << ") = " << i + j << " over square domain [0,1]x[0,1] "
                    << "is not matching using quadrature points.\n";
 
          error_test_2++;
        }
      }
  }
 
  if (n == 1) {
    std::cout << "**********************************\n";
    std::cout << "Triangle Quadrature Test\n";
    std::cout << "**********************************\n";
  }
  std::cout << "Quad order = " << n << ". ";
  if (error_test_1 == 0)
    std::cout << "TEST 1 : PASS. ";
  else
    std::cout << "TEST 1 : FAIL. ";
  if (error_test_2 == 0)
    std::cout << "TEST 2 : PASS. ";
  else
    std::cout << "TEST 2 : FAIL. ";
  std::cout << "\n";
}

References anonymous_namespace{testFeLib.cpp}::checkRefIntegration(), getExactIntegrationRefTri(), anonymous_namespace{testFeLib.cpp}::readElements(), anonymous_namespace{testFeLib.cpp}::readNodes(), anonymous_namespace{testFeLib.cpp}::tol, and util::vtk_type_triangle.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testTriElemTime()

void test::testTriElemTime	(	size_t	n,
		size_t	N
	)

Computes the time needed when quad data for elements are stored and when they are computed as and when needed.

This function allocates dummy elements and test how much time it is required to do computation when the quad data are stored for each element and when the quad data are computed.

Parameters

n	Order of quadrature point approximation
N	Number of elements on which this test is performed

Definition at line 535 of file testFeLib.cpp.

                                             {
 
  // get Quadrature
  auto quad = fe::TriElem(n);
 
  //
  // Test 1
  //
  std::vector<util::Point> nodes = {util::Point(2., 2., 0.),
                                     util::Point(4., 2., 0.),
                                     util::Point(2., 4., 0.)};
  size_t num_vertex = 3;
  //  std::vector<size_t> elements;
  //  for (size_t i = 0; i < 3 * N; i++)
  //    elements.emplace_back(i % 2);
  std::vector<std::vector<size_t>> elements;
  for (size_t i = 0; i < N; i++)
    elements.emplace_back(std::vector<size_t>{0, 1, 2});
 
  auto t11 = steady_clock::now();
  // method 1: Compute quad points at each call
  // loop over elements and compute I_approx
  double sum = 0.;
  for (size_t e = 0; e < N; e++) {
    //    std::vector<util::Point> enodes = {nodes[elements[num_vertex * e +
    //    0]],
    //                                        nodes[elements[num_vertex * e +
    //                                        1]], nodes[elements[num_vertex * e
    //                                        + 2]]};
    std::vector<util::Point> enodes = {
        nodes[elements[e][0]], nodes[elements[e][1]], nodes[elements[e][2]]};
    std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
    for (auto qd : qds)
      sum += qd.d_w * (qd.d_shapes[0] + qd.d_shapes[1] + qd.d_shapes[2]);
  }
  auto t12 = steady_clock::now();
 
  // method 2: Compute quad points in the beginning and use it when needed
  size_t num_quad_pts = 0;
  std::vector<fe::QuadData> quad_data;
  for (size_t e = 0; e < N; e++) {
    std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
    if (e == 0)
      num_quad_pts = qds.size();
    for (auto qd : qds)
      quad_data.emplace_back(qd);
  }
 
  auto t21 = steady_clock::now();
  sum = 0.;
  for (size_t e = 0; e < N; e++) {
    for (size_t q = 0; q < num_quad_pts; q++) {
      fe::QuadData qd = quad_data[e * num_quad_pts + q];
      sum += qd.d_w * (qd.d_shapes[0] + qd.d_shapes[1] + qd.d_shapes[2]);
    }
  }
  auto t22 = steady_clock::now();
 
  if (n == 1 and N == 1000) {
    std::cout << "**********************************\n";
    std::cout << "Quadrature Time Efficiency Test\n";
    std::cout << "**********************************\n";
  }
  std::cout << "Quad order = " << n << ". Num Elements =  " << N << ".\n ";
  double dt_1 = util::methods::timeDiff(t12, t11, "seconds");
  double dt_2 = util::methods::timeDiff(t21, t22, "seconds");
  double perc = (dt_1 - dt_2) * 100. / dt_2;
  double qpt_mem = 13 * sizeof(double);
  double mem2 = double(quad_data.capacity() * qpt_mem) / double(1000000);
  std::cout << "  dt1 = " << dt_1 << ", dt2 = " << dt_2 << ", perc = " << perc
            << ". Mem saved = " << mem2 << " MB.\n";
}

References fe::QuadData::d_shapes, fe::QuadData::d_w, and util::methods::timeDiff().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

testUtilMethods()

void test::testUtilMethods

(

)

Test methods

Definition at line 29 of file testUtilLib.cpp.

                           {
 
  const double tol = 1.e-10;
 
  //
  {
    std::pair<util::Point, util::Point> box = {util::Point(),
                                               util::Point(1., 1., 1.)};
    auto corner_pts = geom::getCornerPoints(3, box);
    auto edges = geom::getEdges(3, box);
    auto xc = geom::getCenter(3, box);
 
    for (size_t i=0; i<2; i++)
      for (size_t j=0; j<2; j++)
        for (size_t k=0; k<2; k++) {
 
          auto p = util::Point(double(i), double(j), double(k));
          bool found_p = false;
          for (auto q : corner_pts) {
            if (q.dist(p) < tol)
              found_p = true;
          }
          if (!found_p)
            errExit(std::format("Error: Can not find corner point {}\n", p.printStr()));
        }
 
    if (xc.dist(util::Point(0.5, 0.5, 0.5)) > tol)
      errExit("Error: getCenter()\n");
  }
 
  //
  {
    if (std::abs(geom::triangleArea(util::Point(0., 0., 0.), util::Point(2., 0., 0.), util::Point(1., 1., 0.)) - 1.) > tol)
      errExit("Error: triangleArea()\n");
  }
 
  //
  {
    std::vector<double> x = {1., 0., 0.};
    std::vector<double> y_check = {1./std::sqrt(2.), -1./std::sqrt(2.), 0.};
    auto y = util::rotateCW2D(x, M_PI * 0.25);
    if (util::methods::l2Dist(y_check, y) > tol)
      errExit("Error: rotateCW2D()\n");
 
    if (util::Point(y_check).dist(util::rotateCW2D(util::Point(x), M_PI * 0.25)) > tol)
      errExit("Error: rotateCW2D()\n");
 
    y_check = {1./std::sqrt(2.), 1./std::sqrt(2.), 0.};
    y = util::rotateACW2D(x, M_PI * 0.25);
    if (util::methods::l2Dist(y_check, y) > tol)
      errExit("Error: rotateACW2D()\n");
  }
 
  //
  {
    auto x = util::Point(1., 0., 0.);
    auto a = util::Point(0., 0., 1.);
    auto y_check = util::Point(0., 1., 0.);
    auto y = util::rotate(x, M_PI * 0.5, a);
    if (y_check.dist(y) > tol)
      errExit(std::format("Error: rotate(). y_check = {}, y = {}\n", y_check.printStr(), y.printStr()));
 
    x = util::Point(1., 1., 1.);
    y_check = util::Point(-1., 1., 1.);
    y = util::rotate(x, M_PI * 0.5, a);
    if (y_check.dist(y) > tol)
      errExit(std::format("Error: rotate(). y_check = {}, y = {}\n", y_check.printStr(), y.printStr()));
  }
 
  //
  {
    auto x1 = util::Point(1., 1., 0.);
    auto x2 = util::Point(1., 0., 0.);
    if (std::abs(M_PI*0.25 - util::angle(x1, x2)) > tol)
      errExit("Error: angle()\n");
 
    x2 = util::Point(0., 0., 1.);
    if (std::abs(M_PI*0.5 - util::angle(x1, x2)) > tol)
      errExit("Error: angle()\n");
 
    x2 = util::Point(0., 1., 1.);
    if (std::abs(M_PI/3. - util::angle(x1, x2)) > tol)
      errExit("Error: angle()\n");
  }
}

References util::angle(), anonymous_namespace{testUtilLib.cpp}::errExit(), geom::getCenter(), geom::getCornerPoints(), geom::getEdges(), util::methods::l2Dist(), util::rotate(), util::rotateACW2D(), util::rotateCW2D(), and geom::triangleArea().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function: