Commit 8e1a4506 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1292 TMP Elasticity case compiles with XTensor only. Many deactivated to make...

#1292 TMP Elasticity case compiles with XTensor only. Many deactivated to make that work though... Seldon has been removed due to conflicts with Blas files also in Xtensor.
parent da1addad
This diff is collapsed.
......@@ -72,6 +72,10 @@
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Elasticity/demo_3d.lua"
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "--overwrite_directory"
isEnabled = "YES">
</CommandLineArgument>
</CommandLineArguments>
</LaunchAction>
<ProfileAction
......
......@@ -40,17 +40,19 @@ namespace MoReFEM
const auto Ndof_per_component = NdofPerComponent();
const auto Nnode = ref_felt.Nnode();
local_node_indexes_ = Wrappers::Seldon::ContiguousSequence(index_first_node_in_elementary_data,
Nnode);
local_dof_indexes_ = Wrappers::Seldon::ContiguousSequence(index_first_dof_in_elementary_data,
Ndof);
local_node_indexes_.resize(Nnode);
std::iota(local_node_indexes_.begin(), local_node_indexes_.end(), index_first_node_in_elementary_data);
local_dof_indexes_.resize(Ndof);
std::iota(local_dof_indexes_.begin(), local_dof_indexes_.end(), index_first_dof_in_elementary_data);
std::vector<int> buf(Ndof_per_component);
for (unsigned int i = 0; i < Ncomponent; ++i)
local_dof_indexes_per_component_.push_back(Wrappers::Seldon
::ContiguousSequence(index_first_dof_in_elementary_data + i * Ndof_per_component,
Ndof_per_component));
{
std::iota(buf.begin(), buf.end(), index_first_dof_in_elementary_data + i * Ndof_per_component);
local_dof_indexes_per_component_.push_back(buf);
}
const auto& ref_felt_space_dim = ref_felt.GetFEltSpaceDimension();
......@@ -76,25 +78,25 @@ namespace MoReFEM
{
auto& ret = ref_felt.GetLocalMatrix<0>(); // mutable...
assert(full_matrix.GetM() >= ret.GetM());
assert(full_matrix.GetN() >= ret.GetN());
assert(full_matrix.shape(0) >= ret.shape(0));
assert(full_matrix.shape(1) >= ret.shape(1));
const auto& row_range = ref_felt.GetLocalNodeIndexList();
const auto ref_felt_space_dim = static_cast<int>(ref_felt.GetFEltSpaceDimension());
#ifndef NDEBUG
Wrappers::Seldon::AssertBlockRowValidity(full_matrix, row_range);
// \todo #1292 Wrappers::Seldon::AssertBlockRowValidity(full_matrix, row_range);
#endif // NDEBUG
assert(static_cast<int>(ref_felt.Nnode()) == row_range.GetM());
assert(ref_felt_space_dim == ret.GetN());
assert(ref_felt.Nnode() == row_range.size());
assert(ref_felt_space_dim == ret.shape(1));
const auto Nrow = row_range.GetM();
const auto Nrow = row_range.size();
for (int i = 0; i < Nrow; ++i)
for (auto i = 0ul; i < Nrow; ++i)
{
auto row_index = row_range(i);
auto row_index = row_range[i];
for (int component = 0; component < ref_felt_space_dim; ++component)
ret(i, component) = full_matrix(row_index, component);
......@@ -109,19 +111,19 @@ namespace MoReFEM
{
auto& ret = ref_felt.GetLocalVector<0>(); // mutable...
assert(full_vector.GetM() >= ret.GetM());
assert(full_vector.shape(0) >= ret.shape(0));
const auto& row_range = ref_felt.GetLocalNodeIndexList();
#ifndef NDEBUG
Wrappers::Seldon::AssertBlockValidity(full_vector, row_range);
// \todo #1292 Wrappers::Seldon::AssertBlockValidity(full_vector, row_range);
#endif // NDEBUG
const auto Nrow = row_range.GetM();
const auto Nrow = row_range.size();
for (int i = 0; i < Nrow; ++i)
for (auto i = 0ul; i < Nrow; ++i)
{
auto row_index = row_range(i);
auto row_index = row_range[i];
ret(i) = full_vector(row_index);
}
......
......@@ -18,7 +18,7 @@
# include <memory>
# include <vector>
# include "ThirdParty/IncludeWithoutWarning/Seldon/Seldon.hpp"
// \todo #1292 # include "ThirdParty/IncludeWithoutWarning/Seldon/Seldon.hpp"
# include "Utilities/LinearAlgebra/Storage/Local/LocalMatrixStorage.hpp"
# include "Utilities/LinearAlgebra/Storage/Local/LocalVectorStorage.hpp"
......@@ -161,7 +161,7 @@ namespace MoReFEM
*
* \return Reference to the list of local node indexes.
*/
const Seldon::Vector<int>& GetLocalNodeIndexList() const noexcept;
const std::vector<int>& GetLocalNodeIndexList() const noexcept;
/*!
* \brief Return the (contiguous) list of dof indexes related to RefFEltInFEltSpace in ElementaryData matrices
......@@ -174,7 +174,7 @@ namespace MoReFEM
*
* \return Reference to the list of local node indexes associated to \a component_index -th component.
*/
const Seldon::Vector<int>& GetLocalDofIndexList(unsigned int component_index) const noexcept;
const std::vector<int>& GetLocalDofIndexList(unsigned int component_index) const noexcept;
/*!
......@@ -185,7 +185,7 @@ namespace MoReFEM
*
* \return Reference to the list of local dof indexes.
*/
const Seldon::Vector<int>& GetLocalDofIndexList() const noexcept;
const std::vector<int>& GetLocalDofIndexList() const noexcept;
//! Returns the number of components.
unsigned int Ncomponent() const noexcept;
......@@ -231,16 +231,16 @@ namespace MoReFEM
* \brief (Contiguous) list of node indexes related to RefFEltInFEltSpace in ElementaryData matrices
* or vectors; this method returns only a subset related to a given component of the unknown.
*/
std::vector<Seldon::Vector<int>> local_dof_indexes_per_component_;
std::vector<std::vector<int>> local_dof_indexes_per_component_;
//! (Contiguous) list of dof indexes related to RefFEltInFEltSpace in ElementaryData matrices or vectors.
Seldon::Vector<int> local_dof_indexes_;
std::vector<int> local_dof_indexes_;
/*!
* \brief Return the (contiguous) list of node indexes related to RefFEltInFEltSpace in ElementaryData matrices
* or vectors; this method returns only a subset related to a given component of the unknown.
*/
Seldon::Vector<int> local_node_indexes_;
std::vector<int> local_node_indexes_;
};
......
......@@ -42,19 +42,19 @@ namespace MoReFEM
}
inline const Seldon::Vector<int>& RefFEltInLocalOperator::GetLocalNodeIndexList() const noexcept
inline const std::vector<int>& RefFEltInLocalOperator::GetLocalNodeIndexList() const noexcept
{
return local_node_indexes_;
}
inline const Seldon::Vector<int>& RefFEltInLocalOperator::GetLocalDofIndexList() const noexcept
inline const std::vector<int>& RefFEltInLocalOperator::GetLocalDofIndexList() const noexcept
{
return local_dof_indexes_;
}
inline const Seldon::Vector<int>& RefFEltInLocalOperator::GetLocalDofIndexList(unsigned int component_index) const noexcept
inline const std::vector<int>& RefFEltInLocalOperator::GetLocalDofIndexList(unsigned int component_index) const noexcept
{
assert(component_index < local_dof_indexes_per_component_.size());
return local_dof_indexes_per_component_[component_index];
......@@ -81,7 +81,7 @@ namespace MoReFEM
inline unsigned int RefFEltInLocalOperator::GetIndexFirstDofInElementaryData() const noexcept
{
assert(static_cast<int>(index_first_dof_) == GetLocalDofIndexList()(0));
assert(static_cast<int>(index_first_dof_) == GetLocalDofIndexList()[0]);
return index_first_dof_;
}
......
......@@ -146,7 +146,7 @@ namespace MoReFEM
scalar_initial_condition_ptr scalar_initial_condition_z_;
//! Content of the parameter.
mutable LocalVector content_;
mutable XtensorVector content_;
};
......
......@@ -68,7 +68,7 @@ namespace MoReFEM
}
inline const LocalVector& ThreeDimensionalInitialCondition::SupplGetValue(const SpatialPoint& coords) const
inline const XtensorVector& ThreeDimensionalInitialCondition::SupplGetValue(const SpatialPoint& coords) const
{
content_(0) = GetScalarInitialConditionX().GetValue(coords);
content_(1) = GetScalarInitialConditionY().GetValue(coords);
......
......@@ -675,8 +675,8 @@ namespace MoReFEM
using initial_condition_type = InputDataNS::InitialCondition<IndexT>;
auto initial_condition_ptr =
Internal::FormulationSolverNS::InitThreeDimensionalInitialCondition<initial_condition_type>(mesh,
input_data);
Internal::FormulationSolverNS::InitThreeDimensionalInitialCondition<initial_condition_type>(mesh,
input_data);
// Nullptr sign the case initial condition should be ignored.
if (!initial_condition_ptr)
......
......@@ -31,13 +31,13 @@ namespace MoReFEM
ComputeJacobian::ComputeJacobian(const unsigned int mesh_dimension)
{
jacobian_.Resize(static_cast<int>(mesh_dimension), static_cast<int>(mesh_dimension));
jacobian_.resize({ mesh_dimension, mesh_dimension });
first_derivate_shape_function_.resize(static_cast<std::size_t>(mesh_dimension));
}
const LocalMatrix& ComputeJacobian::Compute(const GeometricElt& geometric_elt,
const LocalCoords& local_coords)
const XtensorMatrix& ComputeJacobian::Compute(const GeometricElt& geometric_elt,
const LocalCoords& local_coords)
{
const unsigned int Ncoords = geometric_elt.Ncoords();
const unsigned int Ncomponent = geometric_elt.GetDimension();
......@@ -45,13 +45,13 @@ namespace MoReFEM
auto& jacobian = GetNonCstJacobian();
auto& first_derivate_shape_function = GetNonCstFirstDerivateShapeFunction();
const unsigned int Nshape_function = static_cast<unsigned int>(jacobian.GetN());
const unsigned int Nshape_function = static_cast<unsigned int>(jacobian.shape(1));
for (auto& item : first_derivate_shape_function)
item = 0.;
assert(Ncomponent <= Nshape_function);
jacobian.Zero();
jacobian.fill(0.);
for (unsigned int local_coord_index = 0u; local_coord_index < Ncoords; ++local_coord_index)
{
......@@ -83,7 +83,8 @@ namespace MoReFEM
coords_in_geom_elt[component] * first_derivate_shape_function[shape_fct_index];
}
}
// #446
// \todo #446
if (geometric_elt.GetIdentifier() == Advanced::GeometricEltEnum::Point1)
{
jacobian(0, 0) = 1;
......
......@@ -112,8 +112,8 @@ namespace MoReFEM
*
* \return Jacobian matrix.
*/
const LocalMatrix& Compute(const GeometricElt& geom_elt,
const LocalCoords& local_coords);
const XtensorMatrix& Compute(const GeometricElt& geom_elt,
const LocalCoords& local_coords);
private:
......@@ -121,7 +121,7 @@ namespace MoReFEM
* \brief Non constant accessor to the matrix that holds the result of the computation; it should not
* be accessible publicly.
*/
LocalMatrix& GetNonCstJacobian() noexcept;
XtensorMatrix& GetNonCstJacobian() noexcept;
/*!
* \brief Non constant accessor to the values of the first derivates of the shape function,
......@@ -137,7 +137,7 @@ namespace MoReFEM
private:
//! Matrix that holds the result of the computation; it should not be accessible publicly.
LocalMatrix jacobian_;
XtensorMatrix jacobian_;
/*!
* \brief Values of the first derivates of the shape function, updated along \a jacobian_.
......
......@@ -28,7 +28,7 @@ namespace MoReFEM
{
inline LocalMatrix& ComputeJacobian::GetNonCstJacobian() noexcept
inline XtensorMatrix& ComputeJacobian::GetNonCstJacobian() noexcept
{
return jacobian_;
}
......
......@@ -90,76 +90,76 @@ namespace MoReFEM
}
}
LocalCoords Global2local(const GeometricElt& geometric_elt,
const Coords& coords,
const unsigned int mesh_dimension)
{
static_cast<void>(mesh_dimension);
//General Newton routine for the nonlinear case
assert(geometric_elt.GetDimension() == mesh_dimension
&& "This routine works only if the mesh as the same dimension as the elements.");
const unsigned int dimension = geometric_elt.GetDimension();
//Initial point of the routine: barycenter of the reference element.
const auto& ref_geom_elt = geometric_elt.GetRefGeomElt();
LocalCoords current_guess = ref_geom_elt.GetBarycenter();
unsigned int iter = 0u;
unsigned int Niter_max = 100u;
double error = std::numeric_limits<double>::max();
double det;
Coords::unique_ptr buf_ptr = Internal::CoordsNS::Factory::Origin();
Coords::unique_ptr coords_from_guess_ptr = Internal::CoordsNS::Factory::Origin();
auto& buf = *buf_ptr;
auto& coords_from_guess = *coords_from_guess_ptr;
ComputeJacobian compute_jacobian_helper(mesh_dimension);
while ((!NumericNS::IsZero(error)) && (++iter < Niter_max))
{
Local2Global(geometric_elt, current_guess, buf);
decltype(auto) jacobian = compute_jacobian_helper.Compute(geometric_elt, current_guess);
LocalMatrix inv_jacobian = Wrappers::Seldon::ComputeInverseSquareMatrix(jacobian, det);
error = 0.;
for (unsigned int i = 0u; i < dimension; ++i)
{
const int int_i = static_cast<int>(i);
double tmp = 0.;
for (unsigned int j = 0u; j < dimension; ++j)
tmp += inv_jacobian(int_i, static_cast<int>(j)) * (coords[j] - buf[j]);
error += NumericNS::Square(tmp);
current_guess.GetNonCstValue(i) += tmp;
}
error = std::sqrt(error);
}
if (iter == Niter_max)
throw ExceptionNS::GeometricElt::Global2LocalNoConvergence(__FILE__, __LINE__);
{
Local2Global(geometric_elt, current_guess, coords_from_guess);
if (!NumericNS::IsZero(Distance(coords, coords_from_guess)))
throw Exception("Convergence toward a bad solution!", __FILE__, __LINE__);
}
return current_guess;
}
// \todo #1292
// LocalCoords Global2local(const GeometricElt& geometric_elt,
// const Coords& coords,
// const unsigned int mesh_dimension)
// {
// static_cast<void>(mesh_dimension);
//
// //General Newton routine for the nonlinear case
// assert(geometric_elt.GetDimension() == mesh_dimension
// && "This routine works only if the mesh as the same dimension as the elements.");
//
// const unsigned int dimension = geometric_elt.GetDimension();
//
// //Initial point of the routine: barycenter of the reference element.
// const auto& ref_geom_elt = geometric_elt.GetRefGeomElt();
//
// LocalCoords current_guess = ref_geom_elt.GetBarycenter();
//
// unsigned int iter = 0u;
// unsigned int Niter_max = 100u;
// double error = std::numeric_limits<double>::max();
//
// double det;
//
// Coords::unique_ptr buf_ptr = Internal::CoordsNS::Factory::Origin();
// Coords::unique_ptr coords_from_guess_ptr = Internal::CoordsNS::Factory::Origin();
//
// auto& buf = *buf_ptr;
// auto& coords_from_guess = *coords_from_guess_ptr;
//
// ComputeJacobian compute_jacobian_helper(mesh_dimension);
//
// while ((!NumericNS::IsZero(error)) && (++iter < Niter_max))
// {
// Local2Global(geometric_elt, current_guess, buf);
//
// decltype(auto) jacobian = compute_jacobian_helper.Compute(geometric_elt, current_guess);
//
// LocalMatrix inv_jacobian = Wrappers::Seldon::ComputeInverseSquareMatrix(jacobian, det);
//
// error = 0.;
//
// for (unsigned int i = 0u; i < dimension; ++i)
// {
// const int int_i = static_cast<int>(i);
// double tmp = 0.;
//
// for (unsigned int j = 0u; j < dimension; ++j)
// tmp += inv_jacobian(int_i, static_cast<int>(j)) * (coords[j] - buf[j]);
//
// error += NumericNS::Square(tmp);
// current_guess.GetNonCstValue(i) += tmp;
// }
//
// error = std::sqrt(error);
// }
//
// if (iter == Niter_max)
// throw ExceptionNS::GeometricElt::Global2LocalNoConvergence(__FILE__, __LINE__);
//
// {
// Local2Global(geometric_elt, current_guess, coords_from_guess);
//
// if (!NumericNS::IsZero(Distance(coords, coords_from_guess)))
// throw Exception("Convergence toward a bad solution!", __FILE__, __LINE__);
// }
//
//
// return current_guess;
// }
} // namespace GeomEltNS
......
......@@ -30,11 +30,12 @@ namespace MoReFEM
{
const unsigned int geom_elts_size = static_cast<unsigned int>(geom_elts.size());
pseudo_normal_ = std::make_unique<LocalVector>(3);
pseudo_normal_ = std::make_unique<LocalVector>();
pseudo_normal_->resize({ 3 });
auto& pseudo_normal = GetNonCstPseudoNormal();
pseudo_normal.Zero();
pseudo_normal.fill(0.);
unsigned int counter_faces = 0;
......@@ -57,7 +58,7 @@ namespace MoReFEM
if (pseudo_normal_face_ptr != nullptr)
{
Seldon::Add(1., *pseudo_normal_face_ptr, pseudo_normal);
pseudo_normal += *pseudo_normal_face_ptr;
++counter_faces;
}
}
......
......@@ -28,7 +28,8 @@ namespace MoReFEM
{
if (pseudo_normal_ == nullptr)
{
pseudo_normal_ = std::make_unique<LocalVector>(3);
pseudo_normal_ = std::make_unique<LocalVector>();
pseudo_normal_->resize({ 3 });
const unsigned int geom_elts_size = static_cast<unsigned int>(geom_elts.size());
......
......@@ -38,13 +38,14 @@ namespace MoReFEM
{
if (pseudo_normal_ == nullptr)
{
pseudo_normal_ = std::make_unique<LocalVector>(3);
pseudo_normal_ = std::make_unique<LocalVector>();
pseudo_normal_->resize({ 3 });
const unsigned int geom_elts_size = static_cast<unsigned int>(geom_elts.size());
auto& pseudo_normal = GetNonCstPseudoNormal();
pseudo_normal.Zero();
pseudo_normal.fill(0.);
unsigned int index_p1 = 0;
unsigned int index_p2 = 0;
......@@ -113,8 +114,8 @@ namespace MoReFEM
"This should not happen if the mesh is correct.", __FILE__, __LINE__);
const double alpha = std::acos((x * u + y * v + z * w) / (norm_31 * norm_21));
Seldon::Add(alpha, *pseudo_normal_ptr, pseudo_normal);
pseudo_normal += alpha * *pseudo_normal_ptr;
}
}
}
......
......@@ -140,10 +140,10 @@ namespace MoReFEM
const std::vector<unsigned int>& GetVertexCoordsIndexList() const;
//! Constant accessor on the pseudo-normal.
const LocalVector& GetPseudoNormal() const noexcept;
const XtensorVector& GetPseudoNormal() const noexcept;
//! Constant accessor on the pointer of the pseudo-normal.
const std::unique_ptr<LocalVector>& GetPseudoNormalPtr() const noexcept;
const std::unique_ptr<XtensorVector>& GetPseudoNormalPtr() const noexcept;
protected:
......@@ -172,7 +172,7 @@ namespace MoReFEM
void SetVertexCoordsList(const Coords::vector_raw_ptr& coords_list);
//! Non constant accessor on the pseudo-normal.
LocalVector& GetNonCstPseudoNormal() noexcept;
XtensorVector& GetNonCstPseudoNormal() noexcept;
private:
......@@ -212,7 +212,7 @@ namespace MoReFEM
protected:
//! Pseudo-normal of the interface.
std::unique_ptr<LocalVector> pseudo_normal_ = nullptr;
std::unique_ptr<XtensorVector> pseudo_normal_ = nullptr;
};
......
......@@ -93,19 +93,19 @@ namespace MoReFEM
}
inline const LocalVector& Interface::GetPseudoNormal() const noexcept
inline const XtensorVector& Interface::GetPseudoNormal() const noexcept
{
assert(!(!pseudo_normal_));
return *pseudo_normal_;
}
inline LocalVector& Interface::GetNonCstPseudoNormal() noexcept
inline XtensorVector& Interface::GetNonCstPseudoNormal() noexcept
{
return const_cast<LocalVector&>(GetPseudoNormal());
return const_cast<XtensorVector&>(GetPseudoNormal());
}
inline const std::unique_ptr<LocalVector>& Interface::GetPseudoNormalPtr() const noexcept
inline const std::unique_ptr<XtensorVector>& Interface::GetPseudoNormalPtr() const noexcept
{
return pseudo_normal_;
}
......
......@@ -114,7 +114,7 @@ namespace MoReFEM
const auto& pseudo_normal_triangle_ptr = face.GetUnorientedInterface().GetPseudoNormalPtr();
const auto& pseudo_normal_triangle = *pseudo_normal_triangle_ptr;
assert(pseudo_normal_triangle.GetLength() == 3);
assert(pseudo_normal_triangle.size() == 3);
const auto& edge_list = geom_elt.GetOrientedEdgeList();
const auto& vertex_list = geom_elt.GetVertexList();