diff --git a/Unstructured_mesh/Unstructured_mesh.pro.user b/Unstructured_mesh/Unstructured_mesh.pro.user
index c48fb92ecbf550ddde7931730408a09105cc3566..bead3c4b78c2c88405e4c81bde7f3abccd5fdf36 100644
--- a/Unstructured_mesh/Unstructured_mesh.pro.user
+++ b/Unstructured_mesh/Unstructured_mesh.pro.user
@@ -1,6 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE QtCreatorProject>
-<!-- Written by QtCreator 4.10.0, 2019-09-19T11:18:12. -->
+<!-- Written by QtCreator 4.10.0, 2019-09-19T22:58:18. -->
<qtcreator>
<data>
<variable>EnvironmentId</variable>
diff --git a/Unstructured_mesh/main.cpp b/Unstructured_mesh/main.cpp
index f29b9880a7f7356a2f3942988472fa0d4b595b6b..bc90c65cfc0f62fc40c7f23a256c07713b4ae505 100644
--- a/Unstructured_mesh/main.cpp
+++ b/Unstructured_mesh/main.cpp
@@ -164,38 +164,95 @@ void twoPrisms(UnstructuredMesh<3, size_t, double, 6>& mesh3){
DBGCHECK;
}
-template <typename Type, Type startIndex, Type EndIndex, int increment = 1, Type... t>
-struct MakeCustomIntegerSequence : public MakeCustomIntegerSequence<Type, startIndex + increment, EndIndex, increment, t..., startIndex> {
-};
-template <typename Type, Type EndIndex, int increment, Type... t>
-struct MakeCustomIntegerSequence<Type, EndIndex, EndIndex, increment, t...> {
- using type = std::integer_sequence<Type, t..., EndIndex>;
-};
-
-template<typename Type, Type startIndex, Type EndIndex, int increment = 1>
-using make_custom_integer_sequence = typename MakeCustomIntegerSequence<Type, startIndex, EndIndex, increment>::type;
-
-
-template <unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve, unsigned int ... Dimensions>
-MeshDataContainer<Vertex<Dimension, double>, Dimensions...>
-___ComputeCenters(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh, std::integer_sequence<unsigned int, Dimensions...>){
-
- MeshDataContainer<Vertex<Dimension, double>, Dimensions...> centers(mesh);
-
- _ComputeCenters<1, Dimension, Dimensions...>::compute(centers, mesh);
-
- return centers;
-}
+void twoDeformedPrisms(UnstructuredMesh<3, size_t, double, 6>& mesh3){
+ DBGCHECK;
+ mesh3.GetVertices().push_back({0, {0,0,1}});
+ mesh3.GetVertices().push_back({1, {1,0,0}});
+ mesh3.GetVertices().push_back({2, {0,1,0}});
+ mesh3.GetVertices().push_back({3, {1,1,0}});
+ mesh3.GetVertices().push_back({4, {0,0,2}});
+ mesh3.GetVertices().push_back({5, {1,0,1}});
+ mesh3.GetVertices().push_back({6, {0,1,1}});
+ mesh3.GetVertices().push_back({7, {1,1,2}});
+ DBGCHECK;
+ mesh3.GetEdges().push_back({0,0,1});
+ mesh3.GetEdges().push_back({1,0,2});
+ mesh3.GetEdges().push_back({2,1,3});
+ mesh3.GetEdges().push_back({3,2,3});
+ mesh3.GetEdges().push_back({4,0,4});
+ mesh3.GetEdges().push_back({5,1,5});
+ mesh3.GetEdges().push_back({6,3,7});
+ mesh3.GetEdges().push_back({7,2,6});
+ mesh3.GetEdges().push_back({8,4,5});
+ mesh3.GetEdges().push_back({9,5,7});
+ mesh3.GetEdges().push_back({10,6,7});
+ mesh3.GetEdges().push_back({11,4,6});
+ mesh3.GetEdges().push_back({12,6,5});
+ mesh3.GetEdges().push_back({13,2,1});
+ DBGCHECK;
+ size_t index = 0;
+ mesh3.GetFaces().push_back(index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(0,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(1,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(13,true);
+ mesh3.GetFaces().push_back(++index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(13,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(2,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(3,true);
+ mesh3.GetFaces().push_back(++index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(0,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(5,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(8,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(4,true);
+ mesh3.GetFaces().push_back(++index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(1,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(4,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(11,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(7,true);
+ mesh3.GetFaces().push_back(++index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(3,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(6,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(10,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(7,true);
+ mesh3.GetFaces().push_back(++index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(2,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(6,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(9,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(5,true);
+ mesh3.GetFaces().push_back(++index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(8,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(12,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(11,true);
+ mesh3.GetFaces().push_back(++index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(9,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(10,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(12,true);
+ mesh3.GetFaces().push_back(++index);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(12,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(7,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(13,true);
+ mesh3.GetFaces().at(index).GetSubelements().AddSubelement(5,true);
+ DBGCHECK;
+ mesh3.GetFaces().at(0).SetNextBElem(2,0);
+ mesh3.GetFaces().at(2).SetNextBElem(3,0);
+ mesh3.GetFaces().at(3).SetNextBElem(8,0);
+ mesh3.GetFaces().at(8).SetNextBElem(6,0);
+ mesh3.GetFaces().at(6).SetNextBElem(0,0);
+ mesh3.GetCells().push_back(0);
+ mesh3.GetCells().at(0).SetBoundaryElementIndex(0);
-template <unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve>
-auto ComputeCenters(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
+ mesh3.GetFaces().at(1).SetNextBElem(5,1);
+ mesh3.GetFaces().at(5).SetNextBElem(4,1);
+ mesh3.GetFaces().at(4).SetNextBElem(8,1);
+ mesh3.GetFaces().at(8).SetNextBElem(7,1);
+ mesh3.GetFaces().at(7).SetNextBElem(1,1);
+ mesh3.GetCells().push_back(1);
+ mesh3.GetCells().at(1).SetBoundaryElementIndex(1);
- return ___ComputeCenters(mesh, make_custom_integer_sequence<unsigned int, 1, Dimension>{});
+ DBGCHECK;
}
-
-
void testMesh2D() {
using sit = UnstructuredMesh<2, size_t, double>;
sit mesh;
@@ -296,6 +353,21 @@ void testMesh2D() {
for(sit::Cell& cell : mesh.GetCells()){
DBGVAR(centers.GetDataDim<2>().at(cell.GetIndex()))
}
+
+/*
+ DBGMSG("computing measures");
+
+ auto measures = ComputeMeasures(mesh);
+
+
+ for(double edgeM :measures.GetDataDim<1>()) {
+ DBGVAR(edgeM)
+ }
+
+ for(double cellM :measures.GetDataDim<2>()) {
+ DBGVAR(cellM)
+ }
+*/
}
@@ -350,15 +422,29 @@ void testMesh3D() {
//_ComputeCenters<1,3, 3,2,1>::compute<size_t, double, 6>(centers, mesh3);
auto centers = ComputeCenters(mesh3);
- auto& faceCent = centers.GetDataDim<2>();
- for(auto& center : faceCent) {
- DBGVAR(center)
+ for(auto& face : mesh3.GetFaces()) {
+ face.SetCenter(centers.template GetDataDim<2>().at(face.GetIndex()));
+ DBGVAR(face.GetCenter())
}
DBGMSG("cellCenter");
- for(sit3::Cell& cell : mesh3.GetCells()){
- DBGVAR(centers.GetDataDim<3>().at(cell.GetIndex()))
+ for(auto& cell : mesh3.GetCells()) {
+ cell.SetCenter(centers.template GetDataDim<3>().at(cell.GetIndex()));
+ DBGVAR(cell.GetCenter())
+ }
+
+ DBGMSG("measure computation");
+
+ auto measures = ComputeMeasures(mesh3);
+ for(double edgeM : measures.GetDataDim<1>()) {
+ DBGVAR(edgeM)
+ }
+ for(double faceM : measures.GetDataDim<2>()) {
+ DBGVAR(faceM)
}
+ for(double cellM : measures.GetDataDim<3>()) {
+ DBGVAR(cellM)
+ }
//DBGVAR(centers.template GetDataDim<3>().at(0))
@@ -367,6 +453,63 @@ void testMesh3D() {
+void test3DMeshDeformedPrisms() {
+ UnstructuredMesh<3, size_t, double, 6> mesh3;
+ twoDeformedPrisms(mesh3);
+
+ //_ComputeCenters<1,3, 3,2,1>::compute<size_t, double, 6>(centers, mesh3);
+ auto centers = ComputeCenters(mesh3);
+
+ for(auto& face : mesh3.GetFaces()) {
+ face.SetCenter(centers[face]);
+ DBGVAR(face.GetCenter())
+ }
+ DBGMSG("cellCenter");
+ for(auto& cell : mesh3.GetCells()) {
+ cell.SetCenter(centers[cell]);
+ DBGVAR(cell.GetCenter())
+ }
+
+ DBGMSG("measure computation");
+
+ auto measures = ComputeMeasures(mesh3);
+ for(double edgeM : measures.GetDataDim<1>()) {
+ DBGVAR(edgeM)
+ }
+ for(double faceM : measures.GetDataDim<2>()) {
+ DBGVAR(faceM)
+ }
+
+ for(double cellM : measures.GetDataDim<3>()) {
+ DBGVAR(cellM)
+ }
+}
+
+
+
+
+
+
+
+
+void testMeshDataContainer() {
+ UnstructuredMesh<3, size_t, double, 6> mesh3;
+ twoDeformedPrisms(mesh3);
+
+
+ MeshDataContainer<std::tuple<int, double, char, double>, 3,2,0> container(mesh3);
+
+
+
+ for(auto& c : container.GetDataDim<0>()) {
+ c=42;
+ }
+
+ for(auto& v : mesh3.GetVertices()){
+ DBGVAR(container.at(v))
+ }
+}
+
template <unsigned int ... Is>
class ClassA
{
@@ -379,13 +522,53 @@ class ClassA
};
+
+//template <typename ... t> class ClassB{};
+
+template<typename Tuple, unsigned int ... Is>
+class ClassB
+ {
+ public:
+ ClassB (const std::integer_sequence<unsigned int, Is...>, Tuple t)
+ {std::tuple_element_t<0,Tuple> typ = 0;
+ DBGVAR(sizeof... (Is), std::get<0>(std::array<size_t, sizeof...(Is)>{Is...}), std::get<0>(t), typ) }
+
+ };
+
+
+template <typename ... t> class ClassC{};
+
+template<unsigned int ... Is, typename ...Types>
+class ClassC<std::integer_sequence<unsigned int,Is...>, std::tuple<Types...>>
+ {
+ public:
+ ClassC (const std::integer_sequence<unsigned int, Is...>, std::tuple<Types...>)
+ {std::tuple_element_t<0,std::tuple<Types...>> typ = 0;
+ DBGVAR(sizeof... (Is), std::get<0>(std::array<size_t, sizeof...(Is)>{Is...}), typ) }
+
+ ClassC () {
+ std::tuple_element_t<0,std::tuple<Types...>> typ = 42.15;
+ std::tuple_element_t<1,std::tuple<Types...>> typ2 = 42.15;
+ DBGVAR(sizeof... (Is), std::get<0>(std::array<size_t, sizeof...(Is)>{Is...}), typ, typ2)
+ }
+ };
+
+
+
void testTemplate() {
ClassA n(std::make_integer_sequence<unsigned int, 3>{});
UnstructuredMesh<3,size_t, double,6> mesh3;
//MeshDataContainer<Vertex<3, double>, 0,1,2> centers2(mesh3,std::make_integer_sequence<unsigned int, 3>{}, Vertex<3, double>{});
//ComputeCenters(mesh3);
- ClassA u(make_custom_integer_sequence<unsigned int, 10, 0, -2>{});
+ ClassA p(make_custom_integer_sequence<unsigned int, 10, 0, -2>{});
+ std::tuple<double, char> t{};
+ t={1,2};
+ ClassB u(make_custom_integer_sequence<unsigned int, 2, 0, -2>{}, t);
+
+ ClassC<std::integer_sequence<unsigned int, 2,0>, std::tuple<double, char>> c(make_custom_integer_sequence<unsigned int, 2, 0, -2>{}, std::tuple<double, char>{});
+ ClassC<std::integer_sequence<unsigned int, 2,0>, std::tuple<double, char>> cc;
+ ClassC<std::integer_sequence<unsigned int, 2,0>, decltype(std::make_tuple(1.0, 'a'))> ccc;
}
@@ -394,7 +577,9 @@ void testTemplate() {
int main()
{
- testMesh2D();
- testMesh3D();
+ //testMesh2D();
+ //testMesh3D();
+ //test3DMeshDeformedPrisms();
+ testMeshDataContainer();
//testTemplate();
}
diff --git a/Unstructured_mesh/mesh_element.h b/Unstructured_mesh/mesh_element.h
index 4e498035de20bf94e9b00eef2758ff0e00efee69..51214cf8e1e0fe84f99eec3211ccaaf44806fe47 100644
--- a/Unstructured_mesh/mesh_element.h
+++ b/Unstructured_mesh/mesh_element.h
@@ -131,16 +131,7 @@ public:
:MeshElementBase<IndexType>(index), CellBoundaryConnection<IndexType> () {
Subelements.fill({INVALID_INDEX(IndexType), false});
}
- /*
- Vertex<3, Real> ComputeCenter(const MeshElements<3, IndexType, Real, Reserve>& ref){
- Vertex<3,Real> tmpCenter = {};
- for(unsigned char i = 0; i < numberOfElements; i++){
- tmpCenter += ref.template GetElements<1>().at(subElements[i]).ComputeCenter(ref);
- }
- ComputationalySignificantElement<3, Real>::Center = tmpCenter * (1.0 / numberOfElements);
- }
- */
};
@@ -238,21 +229,7 @@ public:
void SetBoundaryElementIndex(IndexType index){
BoundaryElement = index;
}
-/*
- Vertex<MeshDim, Real> ComputeCenter(MeshElements<MeshDim, IndexType, Real, Reserve>& ref){
- IndexType boundary = this->BoundaryElement;
- size_t tmp_boundary = boundary;
- Vertex<MeshDim, Real> tmpCenter = {};
- int numberOfBoundaries = 0;
- do {
- numberOfBoundaries++;
- tmpCenter += ref.GetFaces().at(tmp_boundary).ComputeCenter(ref);
- tmp_boundary = ref.GetFaces().at(tmp_boundary).GetNextBElem(this->GetIndex());
- } while(boundary != tmp_boundary);
- ComputationalySignificantElement<MeshDim, Real>::Center = tmpCenter * (1.0 / numberOfBoundaries);
- return ComputationalySignificantElement<MeshDim, Real>::Center;
- }
-*/
+
};
diff --git a/Unstructured_mesh/mesh_functions.h b/Unstructured_mesh/mesh_functions.h
index ac54ecd5c6532e82bd759f98d641b7240cd2b449..0a273633949282404400147f419355566803e6dc 100644
--- a/Unstructured_mesh/mesh_functions.h
+++ b/Unstructured_mesh/mesh_functions.h
@@ -91,6 +91,18 @@ public:
return data._DataContainer<DataType,pos>::_data;
}
+ template <unsigned int ElementDim, unsigned int Dimension, typename IndexType, typename Real, unsigned int Reserve>
+ DataType& at(MeshElement<Dimension, ElementDim, IndexType, Real, Reserve>& element) {
+ return GetDataDim<ElementDim>().at(element.GetIndex());
+ }
+
+ template <unsigned int ElementDim, unsigned int Dimension, typename IndexType, typename Real, unsigned int Reserve>
+ DataType& operator[](MeshElement<Dimension, ElementDim, IndexType, Real, Reserve>& element) {
+ return GetDataDim<ElementDim>()[element.GetIndex()];
+ }
+
+
+
template <unsigned int Dimension, typename IndexType, typename Real, unsigned int ...Reserve>
MeshDataContainer(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
Alocator<sizeof... (Dimensions) - 1>::AlocateMemory(*this, mesh);
@@ -102,6 +114,159 @@ public:
Alocator<sizeof... (Dimensions) - 1>::AlocateMemory(*this, mesh);
}
};
+
+
+
+
+
+
+/**
+ * @brief The MeshDataContainer struct
+ *
+ * A struct designed to manage data boud to mesh.
+ * Creates a serie of vectors sized acording to dimension.
+ */
+template <typename ...DataTypes, unsigned int ...Dimensions>
+struct MeshDataContainer<std::tuple<DataTypes...>, Dimensions...>{
+private:
+
+ template<unsigned int dim, unsigned int pos, unsigned int _dim>
+ struct DimensionPos : DimensionPos<dim, pos + 1,std::get<pos + 1>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...})>{};
+
+ template<unsigned int dim, unsigned int pos>
+ struct DimensionPos<dim, pos, dim>{
+ static constexpr unsigned int res(){return pos;}
+ };
+
+public:
+ template<unsigned int dim>
+ static constexpr unsigned int DimIndex(){
+ return DimensionPos<dim, 0, std::get<0>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...})>::res();
+ }
+public:
+
+ template<unsigned int pos>
+ using DataType = typename std::tuple_element<pos,std::tuple<DataTypes...>>::type;
+
+ template<unsigned int _Dim, typename Dummy = void>
+ struct _DataContainer : _DataContainer<_Dim - 1, Dummy>{
+ std::vector<DataType<_Dim>> _data;
+ };
+
+ template<typename Dummy>
+ struct _DataContainer<0, Dummy>{
+ std::vector<DataType<0>> _data;
+ };
+
+ template<unsigned int pos, typename dummy = void>
+ struct Alocator{
+ MeshDataContainer<std::tuple<DataTypes...>, Dimensions...>& parent;
+ template<unsigned int Dimension, typename IndexType, typename Real, unsigned int ...Reserve>
+ static void AlocateMemory(MeshDataContainer<std::tuple<DataTypes...>, Dimensions...>& parent ,MeshElements<Dimension, IndexType, Real, Reserve...>& mesh) {
+ parent.template GetDataPos<pos>().resize(
+ mesh.template GetElements<std::get<pos>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...})>().size());
+ Alocator<pos - 1>::AlocateMemory(parent, mesh);
+ }
+ };
+
+ template<typename dummy>
+ struct Alocator<0, dummy>{
+ template<unsigned int Dimension, typename IndexType, typename Real, unsigned int ...Reserve>
+ static void AlocateMemory(MeshDataContainer<std::tuple<DataTypes...>, Dimensions...>& parent ,MeshElements<Dimension, IndexType, Real, Reserve...>& mesh) {
+
+ parent.template GetDataPos<0>().resize(
+ mesh.template GetElements<std::get<0>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...})>().size());
+
+ }
+ };
+
+
+
+
+ /**
+ * @brief data
+ * A structure containing vectors of specified type
+ * alocated to match the mesh elements.
+ */
+ _DataContainer<sizeof... (Dimensions) - 1> data;
+
+
+public:
+
+ /**
+ * @brief GetDataDim
+ * @return
+ */
+ template<unsigned int dim>
+ std::vector<std::tuple_element_t<DimIndex<dim>(), std::tuple<DataTypes...>>>& GetDataDim(){
+ return data._DataContainer<DimIndex<dim>()>::_data;
+ }
+
+
+ /**
+ * @brief GetDataPos
+ * @return
+ */
+ template<unsigned int pos>
+ std::vector<DataType<pos>>& GetDataPos(){
+ return data._DataContainer<pos>::_data;
+ }
+
+ template <unsigned int ElementDim, unsigned int Dimension, typename IndexType, typename Real, unsigned int Reserve>
+ std::tuple_element_t<DimIndex<ElementDim>(), std::tuple<DataTypes...>>& at(MeshElement<Dimension, ElementDim, IndexType, Real, Reserve>& element) {
+ return GetDataDim<ElementDim>().at(element.GetIndex());
+ }
+
+ template <unsigned int ElementDim, unsigned int Dimension, typename IndexType, typename Real, unsigned int Reserve>
+ std::tuple_element_t<DimIndex<ElementDim>(), std::tuple<DataTypes...>>& operator[](MeshElement<Dimension, ElementDim, IndexType, Real, Reserve>& element) {
+ return GetDataDim<ElementDim>()[element.GetIndex()];
+ }
+
+
+
+ template <unsigned int Dimension, typename IndexType, typename Real, unsigned int ...Reserve>
+ MeshDataContainer(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
+ Alocator<sizeof... (Dimensions) - 1>::AlocateMemory(*this, mesh);
+ }
+
+};
+
+
+
+
+
+
+template <typename Type, Type startIndex, Type EndIndex, int increment = 1, Type... t>
+struct MakeCustomIntegerSequence : public MakeCustomIntegerSequence<Type, startIndex + increment, EndIndex, increment, t..., startIndex> {
+};
+
+template <typename Type, Type EndIndex, int increment, Type... t>
+struct MakeCustomIntegerSequence<Type, EndIndex, EndIndex, increment, t...> {
+ using type = std::integer_sequence<Type, t..., EndIndex>;
+};
+
+template<typename Type, Type startIndex, Type EndIndex, int increment = 1>
+using make_custom_integer_sequence = typename MakeCustomIntegerSequence<Type, startIndex, EndIndex, increment>::type;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -129,7 +294,6 @@ struct _ComputeCenters{
DBGMSG(dim);
_ComputeCenters<dim + 1, Dimension, DataDimensions...>::compute(centers, mesh);
}
-
};
template <unsigned int Dimension, unsigned int... DataDimensions>
@@ -180,18 +344,221 @@ struct _ComputeCenters<1, Dimension, DataDimensions...>{
+template <unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve, unsigned int ... Dimensions>
+MeshDataContainer<Vertex<Dimension, double>, Dimensions...>
+___ComputeCenters(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh, std::integer_sequence<unsigned int, Dimensions...>){
+
+ MeshDataContainer<Vertex<Dimension, double>, Dimensions...> centers(mesh);
+
+ _ComputeCenters<1, Dimension, Dimensions...>::compute(centers, mesh);
+
+ return centers;
+}
+
+
template <unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve>
-auto __ComputeCenters(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
-/*
- MeshDataContainer centers(mesh, std::make_integer_sequence<unsigned int, Dimension + 1>{},Vertex<Dimension, Real>{});
+auto ComputeCenters(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
+
+ return ___ComputeCenters(mesh, make_custom_integer_sequence<unsigned int, 1, Dimension>{});
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+template <unsigned int dim, unsigned int Dimension, unsigned int... DataDimensions>
+struct _ComputeMeasures{
+ template <typename IndexType, typename Real, unsigned int ...Reserve>
+ static void compute(MeshDataContainer< Real, DataDimensions...>&,MeshElements<Dimension, IndexType, Real, Reserve...>&){
+ static_assert (Dimension > 3,"The measure computation of mesh of dimension higher than 3 is not implemented yet.");
+ throw std::runtime_error("The measure computation of mesh of dimension higher than 3 is not implemented yet.");
+ }
+};
+
+
+
+template <unsigned int... DataDimensions>
+struct _ComputeMeasures<3, 3, DataDimensions...>{
+ template <typename IndexType, typename Real, unsigned int ...Reserve>
+ static void compute(MeshDataContainer< Real, DataDimensions...>& measures,MeshElements<3, IndexType, Real, Reserve...>& mesh){
+
+ auto& cellMeasures = measures.template GetDataDim<3>();
+
+ for (typename MeshElements<3, IndexType, Real, Reserve...>::template ElemType<3>::type& cell : mesh.GetCells()) {
+ IndexType tmpFace = cell.GetBoundaryElementIndex();
+ Real measure = Real();
+ Vertex<3,Real>& cellCenter = cell.GetCenter();
+ HTMLDBGCOND(cell.GetIndex() == 0, cellCenter[0], cellCenter[1], cellCenter[2]);
+ do {
+ // select 3 different vertices
+ IndexType vAIndex = mesh.GetEdges().at(mesh.GetFaces().at(tmpFace).GetSubelements()[0].index).GetVertexAIndex();
+ IndexType vBIndex = mesh.GetEdges().at(mesh.GetFaces().at(tmpFace).GetSubelements()[0].index).GetVertexBIndex();
+ IndexType vCIndex = mesh.GetEdges().at(mesh.GetFaces().at(tmpFace).GetSubelements()[1].index).GetVertexAIndex();
+ if(vCIndex == vAIndex || vCIndex == vBIndex) {
+ vCIndex = mesh.GetEdges().at(mesh.GetFaces().at(tmpFace).GetSubelements()[1].index).GetVertexBIndex();
+ }
+
+ Vertex<3,Real>& a = mesh.GetVertices().at(vAIndex);
+ Vertex<3,Real>& b = mesh.GetVertices().at(vBIndex);
+ Vertex<3,Real>& c = mesh.GetVertices().at(vCIndex);
+
+ // preparing quiantities
+ Vertex<3,Real> vAmcC = (a-cellCenter);
+ Vertex<3,Real> vBmA = (b-a);
+ Vertex<3,Real> vCmA = (c-a);
+ Real inv_sqrBmA = 1.0 / vBmA.SumOfSquares();
+ Real inv_sqrCmA = 1.0 / vCmA.SumOfSquares();
+
+ Real denominator = 1.0 / (1.0 - (pow(vCmA*vBmA,2) * inv_sqrBmA * inv_sqrCmA));
+
+
+ Real param_t = -denominator * (((vAmcC*vBmA) * inv_sqrBmA) - (inv_sqrBmA*inv_sqrCmA*(vAmcC * vCmA)*(vCmA*vBmA)));
+ //param_t *= inv_sqrBmA;
+ Real param_s = -denominator * (((vAmcC*vCmA) * inv_sqrCmA) - (inv_sqrBmA*inv_sqrCmA*(vAmcC * vBmA)*(vCmA*vBmA)));
+
+ Real distance = (vAmcC + (vBmA * param_t) + (vCmA * param_s)).NormEukleid();
+
+ Real tmp = distance * measures.template GetDataDim<2>().at(tmpFace);
+ measure += tmp / 3.0;
+
+ tmpFace = mesh.GetFaces().at(tmpFace).GetNextBElem(cell.GetIndex());
+ } while (tmpFace != cell.GetBoundaryElementIndex());
+
+ cellMeasures.at(cell.GetIndex()) = measure;
+ }
+ }
+};
+
+template <unsigned int... DataDimensions>
+struct _ComputeMeasures<2, 2, DataDimensions...>{
+ template <typename IndexType, typename Real, unsigned int ...Reserve>
+ static void compute(MeshDataContainer< Real, DataDimensions...>& measures,MeshElements<2, IndexType, Real, Reserve...>& mesh){
+
+ auto& surfaceMeasures = measures.template GetDataDim<2>();
- _ComputeCenters<1, Dimension>::compute(centers, mesh);
+ for (typename MeshElements<2, IndexType, Real, Reserve...>::template ElemType<2>::type& cell : mesh.GetCells()) {
+ IndexType tmpEdge = cell.GetBoundaryElementIndex();
+ Real measure = Real();
+ Vertex<2,Real>& cellCenter = cell.GetCenter();
+ do {
+ Vertex<2,Real>& a = mesh.GetVertices().at(mesh.GetEdges().at(tmpEdge).GetVertexAIndex());
+ Vertex<2,Real>& b = mesh.GetVertices().at(mesh.GetEdges().at(tmpEdge).GetVertexBIndex());
+ double tmp = (cellCenter[0] - a[0]) * (b[1] - a[1]);
+ tmp -= (cellCenter[1] - a[1]) * (b[0] - a[0]);
+ measure += 0.5 * fabs(tmp);
+
+ tmpEdge = mesh.GetEdges().at(tmpEdge).GetNextBElem(cell.GetIndex());
+ } while (tmpEdge != cell.GetBoundaryElementIndex());
+
+ surfaceMeasures.at(cell.GetIndex()) = measure;
+ }
+ }
+};
+
+
+
+template <unsigned int Dimension,unsigned int... DataDimensions>
+struct _ComputeMeasures<2, Dimension, DataDimensions...>{
+ template <typename IndexType, typename Real, unsigned int ...Reserve>
+ static void compute(MeshDataContainer< Real, DataDimensions...>& measures,MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
+
+ auto& surfaceMeasures = measures.template GetDataDim<2>();
+
+ for (typename MeshElements<Dimension, IndexType, Real, Reserve...>::template ElemType<2>::type& face : mesh.template GetElements<2>()) {
+
+ Real measure = Real();
+ Vertex<Dimension,Real>& faceCenter = face.GetCenter();
+ for(auto sube : face.GetSubelements()){
+
+ Vertex<Dimension,Real>& a = mesh.GetVertices().at(mesh.GetEdges().at(sube.index).GetVertexAIndex());
+ Vertex<Dimension,Real>& b = mesh.GetVertices().at(mesh.GetEdges().at(sube.index).GetVertexBIndex());
+
+ Real distance = Real();
+
+ Real param = -1.0*(((a-faceCenter)*(b-a))/((b-a).SumOfSquares()));
+
+ distance = (a-faceCenter+(b-a)*param).NormEukleid();
+
+ Real tmp = distance * measures.template GetDataDim<1>().at(sube.index);
+ measure += tmp * 0.5;
+ }
+ surfaceMeasures.at(face.GetIndex()) = measure;
+ }
+ _ComputeMeasures<3, Dimension, DataDimensions...>::compute(measures, mesh);
+ }
+};
+
+
+
+
+
+
+
+template <unsigned int Dimension, unsigned int... DataDimensions>
+struct _ComputeMeasures<1, Dimension, DataDimensions...>{
+ template <typename IndexType, typename Real, unsigned int ...Reserve>
+ static void compute(MeshDataContainer< Real, DataDimensions...>& measures,MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
- return centers;*/
- return 0;
+ auto& edgeLengths = measures.template GetDataDim<1>();
+
+ for (auto& edge : mesh.GetEdges()) {
+ edgeLengths.at(edge.GetIndex()) = (mesh.GetVertices().at(edge.GetVertexAIndex()) -
+ mesh.GetVertices().at(edge.GetVertexBIndex())).NormEukleid();
+ }
+
+ _ComputeMeasures<2, Dimension, DataDimensions...>::compute(measures, mesh);
+ }
+};
+
+
+
+
+template <unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve, unsigned int ... Dimensions>
+MeshDataContainer<Real, Dimensions...>
+___ComputeMeasures(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh, std::integer_sequence<unsigned int, Dimensions...>){
+
+ MeshDataContainer<Real, Dimensions...> measures(mesh);
+
+ _ComputeMeasures<1, Dimension, Dimensions...>::compute(measures, mesh);
+
+ return measures;
+}
+
+
+template <unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve>
+auto ComputeMeasures(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
+
+ return ___ComputeMeasures(mesh, make_custom_integer_sequence<unsigned int, 1, Dimension>{});
}
+
+
+
+
+
+
+
+
+
+
+
+
+
template<unsigned int MeshDimension, unsigned int ElementDim,typename IndexType, typename Real, unsigned int ...Reserve>
struct CellsVertices {
static void run(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh, IndexType index){
diff --git a/Unstructured_mesh/vertex.h b/Unstructured_mesh/vertex.h
index adac9e449c7c2328cb80dcd6741a120c07b1485f..2821302297e482c3f51fa9f467699095549dd6ef 100644
--- a/Unstructured_mesh/vertex.h
+++ b/Unstructured_mesh/vertex.h
@@ -47,6 +47,15 @@ public:
Vertex<Dim, Real> operator/(const Real&) const;
+ /**
+ * @brief operator *
+ *
+ * Scalar product of two vectors
+ * @param v
+ * @return
+ */
+ Real operator*(const Vertex<Dim, Real>& v);
+
Vertex<Dim, Real>& operator+=(const Vertex<Dim, Real>&);
Vertex<Dim, Real>& operator-=(const Vertex<Dim, Real>&);
@@ -122,6 +131,15 @@ Vertex<Dim, Real> Vertex<Dim, Real>::operator /(const Real& x) const {
return this->operator*(Real(1.0)/x);
}
+
+template<unsigned int Dim, typename Real>
+Real Vertex<Dim, Real>::operator*(const Vertex<Dim, Real> &v)
+{
+ return inlineScalarProduct<Dim, Real>::computation(Coordinates, v.Coordinates);
+}
+
+
+
// Adds value of coordinates of another Point
template <unsigned int Dim, typename Real>
Vertex<Dim, Real>& Vertex<Dim, Real>::operator +=(const Vertex<Dim, Real>& v){
diff --git a/debug/htmllogger.h b/debug/htmllogger.h
index 55655524e2b50640ad821b32e27666b0c705cd3d..c9d6ba7706af2f15df155c7ce26a4e6919679157 100644
--- a/debug/htmllogger.h
+++ b/debug/htmllogger.h
@@ -30,7 +30,6 @@ public:
HtmlLogger(const char* fileName)
:HtmlLogger(){
logFileName = fileName;
- create(fileName);
}
~HtmlLogger(){