diff --git a/CMakeLists.txt b/CMakeLists.txt
index de2b96d..46c2b1f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,9 +2,16 @@ cmake_minimum_required(VERSION 3.12)
project(Shell VERSION 1.0 LANGUAGES CXX)
# Find and load CMake configuration of packages containing this plugin's dependencies
-find_package(Sofa.Component.Controller REQUIRED)
+find_package(Sofa.Config REQUIRED)
+sofa_find_package(Sofa.Component.Controller REQUIRED)
+sofa_find_package(Sofa.Component.Topology.Container.Dynamic REQUIRED)
sofa_find_package(Sofa.Component.StateContainer REQUIRED)
+sofa_find_package(Sofa.Component.Mapping.Linear REQUIRED)
+sofa_find_package(Sofa.GL REQUIRED)
+set(README_FILE README.md)
+
+option(SOFA-PLUGIN_SHELLS_ADAPTIVITY "Enables shells adaptivity" OFF)
# List all files
set(SHELL_SRC_DIR src/Shell)
@@ -13,32 +20,102 @@ set(HEADER_FILES
${SHELL_SRC_DIR}/controller/MeshChangedEvent.h
${SHELL_SRC_DIR}/controller/MeshInterpolator.h
${SHELL_SRC_DIR}/controller/MeshInterpolator.inl
+ ${SHELL_SRC_DIR}/controller/TriangleSwitchExample.h
+ ${SHELL_SRC_DIR}/controller/TriangleSwitchExample.inl
${SHELL_SRC_DIR}/engine/JoinMeshPoints.h
${SHELL_SRC_DIR}/engine/JoinMeshPoints.inl
+ ${SHELL_SRC_DIR}/engine/FindClosePoints.h
+ ${SHELL_SRC_DIR}/engine/FindClosePoints.inl
+ ${SHELL_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.h
+ ${SHELL_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.inl
+ ${SHELL_SRC_DIR}/forcefield/CstFEMForceField.h
+ ${SHELL_SRC_DIR}/forcefield/CstFEMForceField.inl
${SHELL_SRC_DIR}/forcefield/TriangularBendingFEMForceField.h
${SHELL_SRC_DIR}/forcefield/TriangularBendingFEMForceField.inl
+ ${SHELL_SRC_DIR}/forcefield/TriangularShellForceField.h
+ ${SHELL_SRC_DIR}/forcefield/TriangularShellForceField.inl
+ ${SHELL_SRC_DIR}/mapping/BendingPlateMechanicalMapping.h
+ ${SHELL_SRC_DIR}/mapping/BendingPlateMechanicalMapping.inl
+ ${SHELL_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.h
+ ${SHELL_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.inl
+ ${SHELL_SRC_DIR}/misc/PointProjection.h
+ ${SHELL_SRC_DIR}/misc/PointProjection.inl
+ ${SHELL_SRC_DIR}/shells2/fem/BezierShellInterpolation.h
+ ${SHELL_SRC_DIR}/shells2/fem/BezierShellInterpolation.inl
+ ${SHELL_SRC_DIR}/shells2/fem/BezierShellInterpolationM.h
+ ${SHELL_SRC_DIR}/shells2/fem/BezierShellInterpolationM.inl
+ ${SHELL_SRC_DIR}/shells2/forcefield/BezierShellForceField.h
+ ${SHELL_SRC_DIR}/shells2/forcefield/BezierShellForceField.inl
+ ${SHELL_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.h
+ ${SHELL_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.inl
)
+
set(SOURCE_FILES
${SHELL_SRC_DIR}/initShell.cpp
${SHELL_SRC_DIR}/controller/MeshChangedEvent.cpp
${SHELL_SRC_DIR}/controller/MeshInterpolator.cpp
+ ${SHELL_SRC_DIR}/controller/TriangleSwitchExample.cpp
${SHELL_SRC_DIR}/engine/JoinMeshPoints.cpp
+ ${SHELL_SRC_DIR}/engine/FindClosePoints.cpp
+ ${SHELL_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.cpp
+ ${SHELL_SRC_DIR}/forcefield/CstFEMForceField.cpp
${SHELL_SRC_DIR}/forcefield/TriangularBendingFEMForceField.cpp
+ ${SHELL_SRC_DIR}/forcefield/TriangularShellForceField.cpp
+ ${SHELL_SRC_DIR}/mapping/BendingPlateMechanicalMapping.cpp
+ ${SHELL_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.cpp
+ ${SHELL_SRC_DIR}/misc/PointProjection.cpp
+ ${SHELL_SRC_DIR}/shells2/fem/BezierShellInterpolation.cpp
+ ${SHELL_SRC_DIR}/shells2/fem/BezierShellInterpolationM.cpp
+ ${SHELL_SRC_DIR}/shells2/forcefield/BezierShellForceField.cpp
+ ${SHELL_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.cpp
)
-set(README_FILES
- README.md
-)
-# Create the plugin library.
+if(SOFA-PLUGIN_SHELLS_ADAPTIVITY)
+ set(COMPILER_DEFINE "SOFA_BUILD_SHELLS_ADAPTIVITY")
+
+ list(APPEND HEADER_FILES
+ ${SHELL_SRC_DIR}/controller/AdaptiveCuttingController.h
+ ${SHELL_SRC_DIR}/controller/AdaptiveCuttingController.inl
+ ${SHELL_SRC_DIR}/controller/Test2DAdapter.h
+ ${SHELL_SRC_DIR}/controller/Test2DAdapter.inl
+ ${SHELL_SRC_DIR}/misc/Optimize2DSurface.h
+ ${SHELL_SRC_DIR}/misc/Optimize2DSurface.inl
+ ${SHELL_SRC_DIR}/misc/SurfaceParametrization.h
+ ${SHELL_SRC_DIR}/misc/SurfaceParametrization.inl
+ )
+
+ list(APPEND SOURCE_FILES
+ ${SHELL_SRC_DIR}/controller/AdaptiveCuttingController.cpp
+ ${SHELL_SRC_DIR}/controller/Test2DAdapter.cpp
+ ${SHELL_SRC_DIR}/misc/Optimize2DSurface.cpp
+ ${SHELL_SRC_DIR}/misc/SurfaceParametrization.cpp
+ )
+
+ if(SofaGui_FOUND AND SofaOpenglVisual_FOUND)
+ list(APPEND HEADER_FILES
+ ${SHELL_SRC_DIR}/cutting/AdaptiveCutting.h
+ )
+
+ list(APPEND SOURCE_FILES
+ ${SHELL_SRC_DIR}/cutting/AdaptiveCutting.cpp
+ )
+ endif()
+
+endif()
+
+
+# Create the plugin library
add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${README_FILES})
# Link the plugin library to its dependency(ies).
-target_link_libraries(${PROJECT_NAME} Sofa.Component.Controller Sofa.Component.StateContainer)
+target_link_libraries(${PROJECT_NAME}
+ Sofa.Component.Controller
+ Sofa.Component.Topology.Container.Dynamic
+ Sofa.Component.StateContainer
+ Sofa.Component.Mapping.Linear
+ Sofa.GL
+)
-# Create package Config, Version & Target files.
-# Deploy the headers, resources, scenes & examples.
-# Set the plugin 'relocatable' if built within SOFA.
-# --> see SofaMacros.cmake
sofa_create_package_with_targets(
PACKAGE_NAME ${PROJECT_NAME}
PACKAGE_VERSION ${PROJECT_VERSION}
diff --git a/examples/sofapython3/SceneCochlea/__pycache__/Scene_Cochlea.cpython-311.pyc b/examples/sofapython3/SceneCochlea/__pycache__/Scene_Cochlea.cpython-311.pyc
new file mode 100644
index 0000000..3e07ce2
Binary files /dev/null and b/examples/sofapython3/SceneCochlea/__pycache__/Scene_Cochlea.cpython-311.pyc differ
diff --git a/examples/xml/ShellTest.scn b/examples/xml/ShellTest.scn
index 95a6ac0..9af6fcf 100644
--- a/examples/xml/ShellTest.scn
+++ b/examples/xml/ShellTest.scn
@@ -27,11 +27,12 @@
-
+
-
+
+
diff --git a/src/Shell/controller/AdaptiveCuttingController.cpp b/src/Shell/controller/AdaptiveCuttingController.cpp
new file mode 100644
index 0000000..229465f
--- /dev/null
+++ b/src/Shell/controller/AdaptiveCuttingController.cpp
@@ -0,0 +1,43 @@
+#include
+#include
+#include
+
+
+namespace sofa
+{
+
+namespace component
+{
+
+namespace controller
+{
+
+using namespace sofa::defaulttype;
+
+SOFA_DECL_CLASS(AdaptiveCuttingController)
+
+// Register in the Factory
+int AdaptiveCuttingControllerClass = core::RegisterObject(
+ "Controller that handles the cutting method based on mesh adaptivity.")
+#ifdef SOFA_FLOAT
+.add< AdaptiveCuttingController >(true) // default template
+#else
+.add< AdaptiveCuttingController >(true) // default template
+# ifndef SOFA_DOUBLE
+.add< AdaptiveCuttingController >()
+# endif
+#endif
+;
+
+#ifndef SOFA_FLOAT
+template class SOFA_SHELLS_API AdaptiveCuttingController;
+#endif //SOFA_FLOAT
+#ifndef SOFA_DOUBLE
+template class SOFA_SHELLS_API AdaptiveCuttingController;
+#endif //SOFA_DOUBLE
+
+} // namespace controller
+
+} // namespace component
+
+} // namespace sofa
diff --git a/src/Shell/controller/AdaptiveCuttingController.h b/src/Shell/controller/AdaptiveCuttingController.h
new file mode 100644
index 0000000..b1e19a9
--- /dev/null
+++ b/src/Shell/controller/AdaptiveCuttingController.h
@@ -0,0 +1,169 @@
+#ifndef SOFA_COMPONENT_CONTROLLER_ADAPTIVECUTTINGCONTROLLER_H
+#define SOFA_COMPONENT_CONTROLLER_ADAPTIVECUTTINGCONTROLLER_H
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+
+namespace sofa
+{
+
+namespace component
+{
+
+namespace controller
+{
+
+/// Class to shield the data type
+class CuttingAdapter
+{
+public:
+ virtual void setTrackedPoint(const collision::BodyPicked &picked) = 0;
+ virtual void freeTrackedPoint() = 0;
+ virtual void addCuttingPoint() = 0;
+};
+
+template
+class AdaptiveCuttingController : public Controller, public CuttingAdapter
+{
+public:
+ SOFA_CLASS(SOFA_TEMPLATE(AdaptiveCuttingController,DataTypes),Controller);
+
+ typedef typename DataTypes::Coord Coord;
+ typedef typename DataTypes::VecCoord VecCoord;
+ //typedef typename DataTypes::Deriv Deriv;
+ //typedef typename DataTypes::VecDeriv VecDeriv;
+ typedef typename Coord::value_type Real;
+
+ typedef sofa::type::Vec<2, Real> Vec2;
+ typedef sofa::type::Vec<3, Real> Vec3;
+ //typedef sofa::type::Mat<2,2,Real> Mat22;
+ //typedef sofa::type::Mat<3,3,Real> Mat33;
+ //typedef type::vector VecVec2;
+ //typedef type::vector VecVec3;
+
+
+ typedef sofa::core::topology::BaseMeshTopology::Edge Edge;
+ typedef sofa::core::topology::BaseMeshTopology::EdgesAroundVertex EdgesAroundVertex;
+ typedef sofa::component::topology::TriangleSetTopologyContainer::TriangleID Index;
+ typedef sofa::component::topology::TriangleSetTopologyContainer::Triangle Triangle;
+ typedef sofa::component::topology::TriangleSetTopologyContainer::TrianglesAroundVertex TrianglesAroundVertex;
+ typedef sofa::component::topology::TriangleSetTopologyContainer::TrianglesAroundEdge TrianglesAroundEdge;
+ typedef sofa::component::topology::TriangleSetTopologyContainer::EdgesInTriangle EdgesInTriangle;
+ typedef sofa::type::vector VecIndex;
+
+ enum { InvalidID = sofa::core::topology::Topology::InvalidID };
+
+ /// @brief If geometric functinal drops below this value the attached node
+ /// is dropped.
+ Data m_affinity;
+
+ virtual void init();
+ virtual void reinit();
+
+ virtual std::string getTemplateName() const
+ {
+ return templateName(this);
+ }
+
+ static std::string templateName(const AdaptiveCuttingController* = NULL)
+ {
+ return DataTypes::Name();
+ }
+
+ void onEndAnimationStep(const double dt);
+ //void onKeyPressedEvent(core::objectmodel::KeypressedEvent *key);
+
+ void draw(const core::visual::VisualParams* vparams);
+
+
+ void setTrackedPoint(const collision::BodyPicked &picked);
+ void freeTrackedPoint() {
+ // Detach the point
+ m_pointId = InvalidID;
+ // Stop cutting
+ m_cutPoints = 0;
+ m_cutEdge = InvalidID;
+ }
+ void addCuttingPoint();
+
+ /// Whether the cutting is in progress or not.
+ bool cutting() { return m_cutPoints > 0; }
+
+protected:
+
+ AdaptiveCuttingController();
+
+
+private:
+
+ Test2DAdapter* m_adapter;
+ sofa::component::topology::TriangleSetTopologyContainer* m_container;
+ sofa::component::topology::TriangleSetGeometryAlgorithms *m_algoGeom;
+ sofa::component::topology::TriangleSetTopologyAlgorithms *m_algoTopo;
+ sofa::core::behavior::MechanicalState* m_state;
+
+ // TODO: This should go to cutting config (maybe?)
+ bool autoCutting;
+
+
+ /// Closest point in the mstate.
+ Index m_pointId;
+ /// A point on a surface to attract to (valid only if m_pointId != InvalidID).
+ Vec3 m_point;
+ /// Position of m_point projected into rest shape.
+ Vec3 m_pointRest;
+ /// @brief Triangle ID inside which m_point is located (valid only if
+ ///m_pointId != InvalidID).
+ Index m_pointTriId;
+ /// @brief Number of iterations during which the attached node will not be
+ /// reattached.
+ unsigned int m_gracePeriod;
+
+ /// @brief Stored index of the first edge to cut when the first cut has
+ /// been delayed.
+ Index m_cutEdge;
+ /// Last cutting point.
+ Index m_cutLastPoint;
+ /// Cutting operation to perform in this step.
+ VecIndex m_cutList;
+ /// Number of cut points defined.
+ int m_cutPoints;
+
+ void switchPoint(const Vec3 &newPoint, const Index newPointTri,
+ const Index newID, const Index newCutEdge);
+
+ /**
+ * Set edge planned for the next cut.
+ *
+ * @param newCutEdge Index of new edge to use.
+ * @param bKeepProtection Whether to keep edge protection for previous
+ * edge.
+ */
+ void setCutEdge(const Index newCutEdge, const bool bKeepProtection=false);
+
+};
+
+
+} // namespace controller
+
+} // namespace component
+
+} // namespace sofa
+
+#endif // #ifndef SOFA_COMPONENT_CONTROLLER_ADAPTIVECUTTINGCONTROLLER_H
diff --git a/src/Shell/controller/AdaptiveCuttingController.inl b/src/Shell/controller/AdaptiveCuttingController.inl
new file mode 100644
index 0000000..671ad0a
--- /dev/null
+++ b/src/Shell/controller/AdaptiveCuttingController.inl
@@ -0,0 +1,509 @@
+#ifndef SOFA_COMPONENT_CONTROLLER_ADAPTIVECUTTINGCONTROLLER_INL
+#define SOFA_COMPONENT_CONTROLLER_ADAPTIVECUTTINGCONTROLLER_INL
+
+// TODO
+// - protect/unprotect cut edges (m_cutEdge,m_cutList)
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+#define OTHER(x, a, b) ((x == a) ? b : a)
+
+// Return non-zero if triangle with points (a,b,c) is defined in
+// counter-clockwise direction.
+// NOTE: Constrained to 2D!
+#define CCW(a,b,c) (\
+ cross(Vec2(b[0]-a[0], b[1]-a[1]), \
+ Vec2(c[0]-a[0], c[1]-a[1])) > 1e-15)
+
+// Test for intersection between two segments (a,b) and (c,d)
+#define INTERSECT(a,b,c,d) (\
+ (CCW(a,c,d) != CCW(b,c,d)) && (CCW(a,b,c) != CCW(a,b,d)))
+
+namespace sofa
+{
+
+namespace component
+{
+
+namespace controller
+{
+
+template
+AdaptiveCuttingController::AdaptiveCuttingController()
+: m_affinity(initData(&m_affinity, (Real)0.7, "affinity", "Threshold for point attachment (value betwen 0 and 1)."))
+, autoCutting(false)
+, m_pointId(InvalidID)
+, m_pointTriId(InvalidID)
+, m_gracePeriod(0)
+, m_cutEdge(InvalidID)
+, m_cutLastPoint(InvalidID)
+, m_cutPoints(0)
+{
+}
+
+template
+void AdaptiveCuttingController::init()
+{
+ this->getContext()->get(m_adapter);
+ if (m_adapter == NULL) {
+ msg_error() << "Unable to find Test2DAdapter component";
+ return;
+ }
+
+ m_state = dynamic_cast*> (this->getContext()->getMechanicalState());
+ if (!m_state) {
+ msg_error() << "Unable to find MechanicalState";
+ return;
+ }
+
+ this->getContext()->get(m_container);
+ if (m_container == NULL) {
+ msg_error() << "Unable to find triangular topology";
+ return;
+ }
+
+ this->getContext()->get(m_algoGeom);
+ if (m_algoGeom == NULL) {
+ msg_error() << "Unable to find TriangleSetGeometryAlgorithms";
+ return;
+ }
+
+ this->getContext()->get(m_algoTopo);
+ if (m_algoTopo == NULL) {
+ msg_error() << "Unable to find TriangleSetTopologyAlgorithms";
+ return;
+ }
+
+ reinit();
+}
+
+template
+void AdaptiveCuttingController::reinit()
+{
+ *this->f_listening.beginEdit() = true;
+ this->f_listening.endEdit();
+
+ if (m_affinity.getValue() < 0.0 || m_affinity.getValue() > 1.0) {
+ msg_error() << "Affinity must be between 0 and 1.";
+ *m_affinity.beginEdit() = 0.7;
+ m_affinity.endEdit();
+ }
+
+}
+
+
+template
+void AdaptiveCuttingController::onEndAnimationStep(const double /*dt*/)
+{
+ if (m_gracePeriod > 0) m_gracePeriod--;
+
+ const VecCoord& x0 = m_state->read(
+ sofa::core::ConstVecCoordId::restPosition())->getValue();
+
+ // Perform the (delayed) cutting
+ if ((m_cutList.size() > 0) && (m_cutPoints > 2) ) {
+ VecIndex newList, endList;
+ // First remove protection
+ for (VecIndex::const_iterator i=m_cutList.begin();
+ i != m_cutList.end(); i++) {
+ m_adapter->unprotectEdge(*i);
+ }
+ // Then perform the incision
+ bool bReachedBorder;
+ m_algoTopo->InciseAlongEdgeList(m_cutList, newList, endList,
+ bReachedBorder);
+ //m_algoTopo->InciseAlongEdge(m_cutList[0], NULL);
+ m_cutList.clear();
+ }
+
+ // Check the values around attached point and reattach if necessary
+ if ((m_pointId != InvalidID) && !m_gracePeriod && !cutting()) {
+ TrianglesAroundVertex N1 =
+ m_container->getTrianglesAroundVertex(m_pointId);
+ Real min = 1.0;
+ for (unsigned int it=0; itmetricGeom(
+ m_container->getTriangle(N1[it]), N1[it]);
+ if (f < min) min = f;
+ }
+ if (min < m_affinity.getValue()) {
+ // Value too low, try reattaching the node
+ m_gracePeriod = 20;
+ Triangle t = m_container->getTriangle(m_pointTriId);
+ VecIndex pts;
+ for (int i=0; i<3; i++) {
+ if (t[i] != m_pointId) pts.push_back(t[i]);
+ }
+ if ((x0[m_pointId] - x0[ pts[0] ]).norm2() <
+ (x0[m_pointId] - x0[ pts[1] ]).norm2()) {
+ m_pointId = pts[0];
+ } else {
+ m_pointId = pts[1];
+ }
+ }
+ }
+
+
+}
+
+template
+void AdaptiveCuttingController::draw(
+ const core::visual::VisualParams* vparams)
+{
+ //if ((!vparams->displayFlags().getShowBehaviorModels()))
+ return;
+
+ if (!m_state) return;
+ const VecCoord& x = m_state->read(sofa::core::ConstVecCoordId::position())->getValue();
+
+ if (m_cutList.size() > 0) {
+ type::vector points;
+ for (VecIndex::const_iterator i=m_cutList.begin();
+ i != m_cutList.end(); i++) {
+ const Edge &e = m_container->getEdge(*i);
+ points.push_back(x[ e[0] ]);
+ points.push_back(x[ e[1] ]);
+ }
+ vparams->drawTool()->drawLines(points, 4,
+ type::RGBAColor(1.0, 0.0, 0.0, 1.0));
+ }
+
+ if (m_cutEdge != InvalidID) {
+ const Edge &e = m_container->getEdge(m_cutEdge);
+ type::vector points;
+ points.push_back(x[ e[0] ]);
+ points.push_back(x[ e[1] ]);
+ vparams->drawTool()->drawLines(points, 4,
+ type::RGBAColor(1.0, 0.0, 0.0, 1.0));
+ }
+}
+
+
+template
+void AdaptiveCuttingController::setTrackedPoint(
+ const collision::BodyPicked &picked)
+{
+ if (!m_adapter) return;
+
+ using namespace sofa::component::collision;
+
+ //const VecCoord& x0 = m_state->read(
+ // sofa::core::ConstVecCoordId::restPosition())->getValue();
+ const VecCoord& x = m_state->read(
+ sofa::core::ConstVecCoordId::position())->getValue();
+
+
+ // Support only trianglular model! The others don't give any added value.
+ /*if(dynamic_cast(picked.body)) {
+ // Point
+ m_pointId = picked.indexCollisionElement;
+ } else if(dynamic_cast(picked.body)) {
+ // Edge
+ Edge e = m_container->getEdge(picked.indexCollisionElement);
+ Real d1 = (x[ e[0] ] - m_point).norm2();
+ Real d2 = (x[ e[1] ] - m_point).norm2();
+ m_pointId = (d1 < d2 ? e[0] : e[1]);
+ } else*/
+ if(dynamic_cast(picked.body)) {
+
+ Index newId = InvalidID;
+ Index newCutEdge = InvalidID;
+
+ if (!cutting()) {
+ Triangle t = m_container->getTriangle(
+ picked.indexCollisionElement);
+ Real d1 = (x[ t[0] ] - picked.point).norm2();
+ Real d2 = (x[ t[1] ] - picked.point).norm2();
+ Real d3 = (x[ t[2] ] - picked.point).norm2();
+ newId = (d1 < d2) ?
+ (d1 < d3 ? t[0] : t[2]) :
+ (d2 < d3 ? t[1] : t[2]);
+ } else {
+ // During cutting we should allow change only to point directly
+ // connected with last cut point.
+ // We pick a point from N1 shell which is closest to the tracked point.
+
+ // NOTE: We don't allow fixed/boundary nodes (now)
+
+
+ Triangle t = m_container->getTriangle(
+ picked.indexCollisionElement);
+ if ((t[0] == m_cutLastPoint) || (t[1] == m_cutLastPoint) ||
+ (t[2] == m_cutLastPoint)) {
+ // The point is inside a triangle in N1-ring, take one of it's
+ // corner nodes
+ Index pt1 = OTHER(m_cutLastPoint, t[0], t[1]),
+ pt2 = OTHER(m_cutLastPoint, t[2], t[1]);
+
+
+ if (!m_adapter->isPointNormal(pt1)) {
+ if (!m_adapter->isPointNormal(pt2)) {
+ // No point available!
+ newId = InvalidID;
+ } else {
+ newId = pt2;
+ }
+ } else if (!m_adapter->isPointNormal(pt2)) {
+ newId = pt1;
+ } else if (
+ (x[pt1] - picked.point).norm2() <
+ (x[pt2] - picked.point).norm2()) {
+ newId = pt1;
+ } else {
+ newId = pt2;
+ }
+
+ } else {
+
+ // Tracked position is outside the N1-ring. Pick closest point
+ // while checking for crossing segments -- between last cut
+ // edge and (seleted-picked). Idealy we should check with all
+ // cut edges, but that's too much work.
+
+ EdgesAroundVertex N1e =
+ m_container->getEdgesAroundVertex(m_cutLastPoint);
+
+ Index lastCut = InvalidID;
+ Edge ce(InvalidID, InvalidID);
+ if (m_cutList.size() > 0) {
+ lastCut = m_cutList.back();
+ ce = m_container->getEdge(lastCut);
+ }
+
+ Real minDist = DBL_MAX;
+
+ for (Index ie=0; iegetEdge(N1e[ie]);
+ Index otherPt = OTHER(m_cutLastPoint, e[0], e[1]);
+
+ if (!m_adapter->isPointNormal(otherPt))
+ continue;
+
+ Real dist = (picked.point - x[otherPt]).norm2();
+ bool inter = false;
+ if (lastCut != InvalidID) {
+ inter = INTERSECT(x[ce[0]], x[ce[1]],
+ x[otherPt], picked.point);
+ }
+ if ((dist < minDist) && !inter) {
+ minDist = dist;
+ newId = otherPt;
+ newCutEdge = N1e[ie];
+ }
+ }
+
+ }
+ if (newId != InvalidID) {
+ newCutEdge = m_container->getEdgeIndex(m_cutLastPoint, newId);
+ }
+ }
+
+ if (newId == InvalidID) {
+ msg_error() << "Failed to pick a point!";
+ } else {
+ switchPoint(picked.point, picked.indexCollisionElement, newId,
+ newCutEdge);
+ }
+ } else {
+ freeTrackedPoint();
+ }
+
+
+ // Add new cut point during automated cutting
+ if (cutting() && autoCutting) {
+ // Compare edge lengths in triangles sharing the cut edge.
+
+ TrianglesAroundEdge tris =
+ m_container->getTrianglesAroundEdge(m_cutEdge);
+ Real edgeSum = 0.0;
+ int edgeCount = 0;
+
+ for (int i=0; i<2; i++) {
+ EdgesInTriangle elist = m_container->getEdgesInTriangle(tris[i]);
+ for (int j=0; j<3; j++) {
+ if (elist[j] == m_cutEdge) continue;
+ Edge e = m_container->getEdge(elist[j]);
+ edgeSum += (x[e[0]] - x[e[1]]).norm2();
+ edgeCount++;
+ }
+ }
+
+ // Compute mean of (squared) edge lengths
+ if (edgeCount > 0) {
+ edgeSum /= edgeCount;
+
+ if (edgeSum < (x[m_cutLastPoint] - x[m_pointId]).norm2()) {
+ addCuttingPoint();
+ }
+ }
+ }
+}
+
+template
+void AdaptiveCuttingController::addCuttingPoint()
+{
+ if (!m_adapter) return;
+ if (!m_algoTopo || !m_algoGeom) return;
+
+ if (m_pointId == InvalidID) {
+ msg_error() << "BUG! Attempted cutting with no point tracked.";
+ return;
+ }
+
+ const VecCoord& x = m_state->read(
+ sofa::core::ConstVecCoordId::position())->getValue();
+ //const VecCoord& xrest = m_state->read(
+ // sofa::core::ConstVecCoordId::restPosition())->getValue();
+ Coord oldpos = x[m_pointId];
+
+ bool bFirst = !cutting();
+
+ // Check if the target location is in one of the triangles connected to
+ // poin m_pointId.
+ bool bConnected = false;
+ Triangle t = m_container->getTriangle(m_pointTriId);
+ for (int i=0; i<3; i++) {
+ if (t[i] == m_pointId) {
+ bConnected = true;
+ break;
+ }
+ }
+
+ // Make sure the tracked point is at the target location.
+ if (!bConnected) {
+ // NOTE: We can try adding the cut point as far as possible, then
+ // reattach and repeat. But we risk "eating" up all available
+ // points quitckly. Another possibility is waiting a few iterations
+ // before adding another cut point, but this will only work if the
+ // cut point is not last (end of the cut).
+ msg_error() << "Failed to insert cut point! Tracking too slow.";
+ return;
+ } else if (m_adapter->isPointNormal(m_pointId)) {
+ m_adapter->relocatePoint(m_pointId, m_point, m_pointTriId, false);
+ } else {
+ // TODO: We need a parameter controling how close one has to be to the
+ // boundary/fixed node to attach to it. For boundary nodes we have to
+ // project the target position onto the boundary nodes. Fixed points,
+ // of course, have to be kept intact. Or alternatively we may prevent
+ // attachment to fixed nodes completely.
+ msg_error() << "Handling of boundary/fixed nodes is not implemented!";
+ }
+
+
+ // Chose another point in the direction of movement.
+ // TODO: the following is not good
+ Vec3 dir = m_point - oldpos;
+ dir.normalize();
+ // END_TODO
+
+ // Get another point in that direction.
+ Index tId = m_algoGeom->getTriangleInDirection(m_pointId, dir);
+ if (tId == InvalidID) {
+ msg_error() << "BUG! Nothing in cutting direction!";
+ return;
+ }
+
+ EdgesInTriangle elist = m_container->getEdgesInTriangle(tId);
+ Real alpha[2];
+ Index otherPt[2], otherEdge[2];
+ for (int i=0,j=0; i<3; i++) {
+ Edge e = m_container->getEdge(elist[i]);
+ if ((e[0] != m_pointId) && (e[1] != m_pointId)) continue;
+ Vec3 edgeDir;
+ if (e[0] == m_pointId) {
+ otherPt[j] = e[1];
+ edgeDir = x[ e[1] ] - x[ e[0] ];
+ } else {
+ otherPt[j] = e[0];
+ edgeDir = x[ e[0] ] - x[ e[1] ];
+ }
+ otherEdge[j] = elist[i];
+ edgeDir.normalize();
+ alpha[j] = dir * edgeDir;
+ j++;
+ }
+
+ Index newPt = (alpha[0] < alpha[1]) ? otherPt[1] : otherPt[0];
+ Index oldPt = m_pointId;
+ Index edge = (alpha[0] < alpha[1]) ? otherEdge[1] : otherEdge[0];
+
+ m_cutPoints++;
+
+ // Attach the new point and move it to be near the cursor.
+ m_pointId = newPt;
+ m_pointTriId = tId;
+ m_gracePeriod = 20;
+ m_adapter->relocatePoint(newPt,
+ x[oldPt] + (dir*m_adapter->getPrecision()/2.0),
+ tId, false);
+
+ m_cutLastPoint = oldPt;
+ if (!bFirst && (m_cutEdge != InvalidID)) {
+ m_cutList.push_back(m_cutEdge);
+ }
+ setCutEdge(edge, true);
+ // Fix the point so nothing happens to it before the topological change
+ // occurs.
+ m_adapter->setPointFixed(oldPt);
+}
+
+template
+void AdaptiveCuttingController::switchPoint(const Vec3 &newPoint,
+ const Index newPointTri, const Index newId, const Index newCutEdge)
+{
+ if (newId == InvalidID) {
+ freeTrackedPoint();
+ return;
+ }
+
+ if (!m_gracePeriod && (newId != m_pointId)) {
+ m_gracePeriod = 5; // NOTE: Don't put too large value here, or we
+ // will fail to follow quick changes.
+ m_pointId = newId;
+ if (newCutEdge != InvalidID) {
+ setCutEdge(newCutEdge);
+ }
+ }
+ //if (newId == m_pointId) {
+ m_point = newPoint;
+ m_pointTriId = newPointTri;
+ //}
+
+ m_adapter->setPointAttraction(m_pointId, m_point, m_pointTriId);
+
+}
+
+template
+void AdaptiveCuttingController::setCutEdge(const Index newCutEdge,
+ const bool bKeepProtection)
+{
+ if (!m_adapter) return;
+
+ if (!bKeepProtection) {
+ m_adapter->unprotectEdge(m_cutEdge);
+ }
+ m_adapter->protectEdge(newCutEdge);
+
+ m_cutEdge = newCutEdge;
+}
+
+} // namespace controller
+
+} // namespace component
+
+} // namespace sofa
+
+
+#undef OTHER
+#undef CCW
+#undef INTERSECT
+
+
+#endif // #ifndef SOFA_COMPONENT_CONTROLLER_ADAPTIVECUTTINGCONTROLLER_INL
diff --git a/src/Shell/controller/CudaTest2DAdapter.cpp b/src/Shell/controller/CudaTest2DAdapter.cpp
new file mode 100644
index 0000000..26fce42
--- /dev/null
+++ b/src/Shell/controller/CudaTest2DAdapter.cpp
@@ -0,0 +1,46 @@
+
+#include
+#include
+#include
+#include
+
+
+namespace sofa
+{
+
+namespace component
+{
+
+namespace controller
+{
+
+template class Test2DAdapter;
+//#ifdef SOFA_GPU_CUDA_DOUBLE
+//template class Test2DAdapter;
+//#endif // SOFA_GPU_CUDA_DOUBLE
+
+} // namespace controller
+
+} // namespace component
+
+namespace gpu
+{
+
+namespace cuda
+{
+
+SOFA_DECL_CLASS(CudaTest2DAdapter)
+
+// Register in the Factory
+int CudaTest2DAdapterClass = core::RegisterObject("Adaptive mesh improvement component for 2D triangular meshes (for testing) -- the CUDA version")
+.add< component::controller::Test2DAdapter >()
+//#ifdef SOFA_GPU_CUDA_DOUBLE
+//.add< component::controller::Test2DAdapter >()
+//#endif // SOFA_GPU_CUDA_DOUBLE
+;
+
+} // namespace cuda
+
+} // namespace gpu
+
+} // namespace sofa
diff --git a/src/Shell/controller/CudaTest2DAdapter.cu b/src/Shell/controller/CudaTest2DAdapter.cu
new file mode 100644
index 0000000..3852824
--- /dev/null
+++ b/src/Shell/controller/CudaTest2DAdapter.cu
@@ -0,0 +1,613 @@
+#include
+#include
+#include
+//#include
+
+#if defined(__cplusplus) && CUDA_VERSION < 2000
+namespace sofa
+{
+namespace gpu
+{
+namespace cuda
+{
+#endif
+
+extern "C"
+{
+void Test2DAdapterCuda3f_computeTriangleNormal(unsigned int size, const void* x, void* tri);
+void Test2DAdapterCuda3f_functionalGeom(unsigned int size, const void* x, void* tri);
+void Test2DAdapterCuda3f_reduceStep(unsigned int size, void* x, void* pt, void* indices);
+void Test2DAdapterCuda3f_restoreUnchanged(unsigned int size, void* x, void* pt, void* indices);
+void Test2DAdapterCuda3f_smooth(unsigned int size, void* x, void* tri, void* pt, void* indices);
+void Test2DAdapterCuda3f_testAcceptable(unsigned int size, void* x, const void* tri, void* pt, void* indices, float tolerance);
+///// parallel
+void Test2DAdapterCuda3f_prepareGradients(unsigned int size, const void* x, void* tri);
+void Test2DAdapterCuda3f_smoothParallel(unsigned int size, void* x, const void* tri, void* pt);
+void Test2DAdapterCuda3f_reduceStepP(unsigned int size, void* x, void* pt);
+void Test2DAdapterCuda3f_testAcceptableP(unsigned int size, void* x, const void* tri, void* pt, float tolerance);
+void Test2DAdapterCuda3f_restoreUnchangedP(unsigned int size, void* x, void* pt);
+//#ifdef SOFA_GPU_CUDA_DOUBLE
+//void Test2DAdapterCuda3d_computeTriangleNormal(const void* x, const void* n);
+//#endif
+}// extern "C"
+
+typedef unsigned int Index;
+
+// NOTE: must be equivalent to the Test2DAdapterData::TriangleData
+template
+struct TriangleData {
+ Index nodes[3];
+ CudaVec3 normal;
+ Real functional;
+ CudaVec3 gradient[3];
+};
+
+// NOTE: must be equivalent to the Test2DAdapterData::PointData
+template
+struct PointData {
+ bool bBoundary;
+
+ unsigned int nNeighboursPt;
+ const Index *neighboursPt;
+
+ unsigned int nNeighboursTri;
+ const Index *neighboursTri;
+
+ bool bAccepted; /// New position has been accepted in current step.
+ CudaVec3 oldpos;
+ Real oldworst;
+ Real newworst;
+
+ Index mintri;
+ CudaVec3 grad;
+};
+
+//////////////////////
+// GPU-side methods //
+//////////////////////
+
+template
+__device__ CudaVec3 computeTriangleNormal(const CudaVec3* x,
+ const Index nodes[3])
+{
+ CudaVec3 A, B;
+ A = x[ nodes[1] ] - x[ nodes[0] ];
+ B = x[ nodes[2] ] - x[ nodes[0] ];
+
+ CudaVec3 normal = CudaVec3::make(0.0, 0.0, 0.0);
+
+ Real An = invnorm(A), Bn = invnorm(B);
+ if (An > 1e-20 && Bn > 1e-20) {
+ A = A*An;
+ B = B*Bn;
+ normal = cross(A, B);
+ normal = normal * invnorm(normal);
+ }
+
+ return normal;
+}
+
+template
+__device__ Real getMinFunc(Index v, const TriangleData* tri,
+ const PointData* pt)
+{
+ unsigned int nElem = pt[v].nNeighboursTri;
+ Real value = 1.0;
+ // TODO: do some unrolling?
+ for (Index it=0; it tri[ pt[v].neighboursTri[it] ].functional) {
+ value = tri[ pt[v].neighboursTri[it] ].functional;
+ }
+ }
+
+ return value;
+}
+
+template
+__device__ Real functionalGeom(const Index t,
+ const CudaVec3* x, const TriangleData* tri)
+{
+ // TODO: move outside, pass nodes as arguments
+
+ // TODO: we can precompute these, is it worth it?
+ CudaVec3 ab = x[ tri[t].nodes[1] ] - x[ tri[t].nodes[0] ];
+ CudaVec3 ca = x[ tri[t].nodes[0] ] - x[ tri[t].nodes[2] ];
+ CudaVec3 cb = x[ tri[t].nodes[1] ] - x[ tri[t].nodes[2] ];
+
+ // Normalizing factor so that the value is 1 in maximum
+ // TODO: does compiler precompute this? NOTE: is float
+ Real m = 2 * sqrt(3.0f);
+
+ m *= norm(cross(ca,cb)); // || CA × CB ||
+ m /= norm2(ca) + norm2(ab) + norm2(cb);
+
+ // Is triangle inverted?
+ CudaVec3 nnew = computeTriangleNormal(x, tri[t].nodes);
+ if (dot(nnew, tri[t].normal) < 0.0) {
+ m *= -1.0;
+ }
+
+ return m;
+}
+
+template
+__device__ __inline__ Index translateIndexInTriangle(Index index,
+ const TriangleData &tri)
+{
+ if (tri.nodes[0] == index) return 0;
+ if (tri.nodes[1] == index) return 1;
+ return 2;
+}
+
+//////////////////////
+// Kernels //
+//////////////////////
+
+template
+__global__ void Test2DAdapterCuda3t_computeTriangleNormal_kernel(unsigned int size,
+ const CudaVec3* x, TriangleData* tri)
+{
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+ tri[index].normal = computeTriangleNormal(x, tri[index].nodes);
+ }
+}
+
+template
+__global__ void Test2DAdapterCuda3t_functionalGeom_kernel(unsigned int size,
+ const CudaVec3* x, TriangleData* tri)
+{
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+ tri[index].functional = functionalGeom(index, x, tri);
+ }
+}
+
+
+// Laplacian smoothing
+template
+__global__ void Test2DAdapterCuda3t_smoothLaplacian_kernel(unsigned int size,
+ CudaVec3* x, const TriangleData* tri, PointData* pt,
+ const Index *indices)
+{
+ typedef CudaVec3 Vec3;
+
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+
+ Index v = indices[index];
+
+ pt[v].oldpos = x[v];
+ pt[v].oldworst = getMinFunc(v, tri, pt);
+ pt[v].bAccepted = false;
+
+ // Compute centroid of polygon from 1-ring around the vertex
+ Vec3 xnew = Vec3::make(0.0, 0.0, 0.0);
+ for (Index ie=0; ie
+__global__ void Test2DAdapterCuda3t_smoothOptimize_kernel(unsigned int size,
+ CudaVec3* x, const TriangleData* tri, PointData* pt,
+ const Index *indices)
+{
+ typedef CudaVec3 Vec3;
+
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+
+ Index v = indices[index];
+ Vec3 xold = x[v];
+
+ pt[v].oldpos = x[v];
+ pt[v].oldworst = getMinFunc(v, tri, pt);
+ pt[v].bAccepted = false;
+
+ unsigned int nElem = pt[v].nNeighboursTri;
+
+ // Compute gradients
+ // TODO: do it once for all elements!
+ Vec3 grad[10]; // TODO: Vec3 grad[nElem];
+ if (nElem > 10) nElem = 10; // XXX
+ Real delta = 1e-5;
+
+ // NOTE: Constrained to 2D!
+ // TODO: can we use shared memory here?
+ // -- X
+ x[v].x += delta;
+ for (Index it=0; it(pt[v].neighboursTri[it], x, tri);
+ grad[it].x = (m - tri[ pt[v].neighboursTri[it] ].functional)/delta;
+ }
+ // -- Y
+ x[v].x = xold.x;
+ x[v].y += delta;
+ for (Index it=0; it(pt[v].neighboursTri[it], x, tri);
+ grad[it].y = (m - tri[ pt[v].neighboursTri[it] ].functional)/delta;
+ }
+
+ // Find smallest functional with non-zero gradient
+ Index imin = 0;
+ Real fmin = 1.0;
+ for (Index it=0; it 1e-15)) {
+ fmin = tri[ pt[v].neighboursTri[it] ].functional;
+ imin = it;
+ }
+ }
+
+ Vec3 step = grad[imin];
+ // Find out step size
+ Real gamma = 0.05;
+ //gamma *= step.norm();
+ step = step * invnorm(step);
+
+ x[v] = xold + gamma*step;
+ }
+}
+
+template
+__global__ void Test2DAdapterCuda3t_reduceStep_kernel(unsigned int size,
+ CudaVec3* x, const PointData* pt, const Index *indices)
+{
+ typedef CudaVec3 Vec3;
+
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+
+ Index v = indices[index];
+ if (!pt[v].bAccepted) {
+ // The correct step size is best found empiricaly
+ x[v] = pt[v].oldpos + (x[v] - pt[v].oldpos) * Real(2.0/3.0);
+ //x[v] = (x[v] + pt[v].oldpos)/2.0;
+ }
+ }
+}
+
+template
+__global__ void Test2DAdapterCuda3t_restoreUnchanged_kernel(unsigned int size,
+ CudaVec3* x, const PointData* pt, const Index *indices)
+{
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+ Index v = indices[index];
+ if (!pt[v].bAccepted) {
+ x[v] = pt[v].oldpos;
+ }
+ }
+}
+
+template
+__global__ void Test2DAdapterCuda3t_testAcceptable_kernel(unsigned int size,
+ CudaVec3* x, const TriangleData* tri, PointData* pt, const Index *indices, float tolerance)
+{
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+
+ //// This check is not worth the effort
+ //if ((xold - x[v]).norm2() < 1e-8) {
+ // // No change in position
+ // //std::cout << "No change in position for " << v << "\n";
+ // break;
+ //}
+
+ Index v = indices[index];
+ //if (!pt[v].bAccepted) { // TODO
+
+ // We accept any change that doesn't decrease worst metric for the
+ // triangle set.
+ Real newworst = getMinFunc(v, tri, pt);
+ if (newworst >= (pt[v].oldworst + tolerance)) {
+ pt[v].bAccepted = true;
+ }
+ pt[v].newworst = newworst;
+
+ //}
+ }
+}
+
+///// parallel
+
+template
+__global__ void Test2DAdapterCuda3t_prepareGradients_kernel(unsigned int size,
+ CudaVec3* x, TriangleData* tri)
+{
+ typedef CudaVec3 Vec3;
+
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+ // TODO: handle boundary vertices
+
+ for (int i=0; i<3; i++) {
+ // For each corner
+
+ Real delta = 1e-5;
+
+ Index v = tri[index].nodes[i];
+ Vec3 xold = x[v];
+
+ // NOTE: Constrained to 2D!
+ // TODO: can we use shared memory here?
+ // -- X
+ x[v].x += delta;
+ Real m = functionalGeom(index, x, tri);
+ tri[index].gradient[i].x = (m - tri[index].functional)/delta;
+ // -- Y
+ x[v].x = xold.x;
+ x[v].y += delta;
+ m = functionalGeom(index, x, tri);
+ tri[index].gradient[i].y = (m - tri[index].functional)/delta;
+
+ x[v].y = xold.y;
+ }
+ }
+}
+
+template
+__global__ void Test2DAdapterCuda3t_smoothParallel_kernel(unsigned int size,
+ CudaVec3* x, const TriangleData* tri, PointData* pt)
+{
+ typedef CudaVec3 Vec3;
+
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+
+ pt[index].oldpos = x[index];
+ pt[index].oldworst = getMinFunc(index, tri, pt);
+ pt[index].bAccepted = false;
+
+ unsigned int nElem = pt[index].nNeighboursTri;
+
+ // Find smallest functional with non-zero gradient
+ Index tmin = Index(-1);
+ Real fmin = 1.0;
+ Vec3 step = Vec3::make(0.0, 0.0, 0.0);
+ if (!pt[index].bBoundary) {
+ for (Index it=0; it 1e-15)) {
+ fmin = tri[t].functional;
+ tmin = t;
+ step = tg;
+ }
+ }
+ }
+ pt[index].mintri = tmin;
+ pt[index].grad = step;
+
+ // Sync -- we need mintri to be available for all points
+ __syncthreads();
+
+ // Consult neighbourhood and make an estimate
+ Vec3 ns = Vec3::make(0.0, 0.0, 0.0);
+ if (!pt[index].bBoundary) {
+
+ for (Index it=0; it
+__global__ void Test2DAdapterCuda3t_reduceStepP_kernel(unsigned int size,
+ CudaVec3* x, const PointData* pt)
+{
+ typedef CudaVec3 Vec3;
+
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+
+ if (!pt[index].bAccepted) {
+ // The correct step size is best found empiricaly
+ x[index] = pt[index].oldpos
+ + (x[index] - pt[index].oldpos) * Real(2.0/3.0);
+ //x[index] = (x[v] + pt[index].oldpos)/2.0;
+ }
+ }
+}
+
+template
+__global__ void Test2DAdapterCuda3t_testAcceptableP_kernel(unsigned int size,
+ CudaVec3* x, const TriangleData* tri, PointData* pt,
+ float tolerance)
+{
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+
+ //// This check is not worth the effort
+ //if ((xold - x[index]).norm2() < 1e-8) {
+ // // No change in position
+ // //std::cout << "No change in position for " << index << "\n";
+ // break;
+ //}
+
+ //if (!pt[index].bAccepted) { // TODO
+
+ // We accept any change that doesn't decrease worst metric for the
+ // triangle set.
+ Real newworst = getMinFunc(index, tri, pt);
+ if (newworst >= (pt[index].oldworst + tolerance)) {
+ pt[index].bAccepted = true;
+ }
+ pt[index].newworst = newworst;
+
+ //}
+ }
+}
+
+template
+__global__ void Test2DAdapterCuda3t_restoreUnchangedP_kernel(unsigned int size,
+ CudaVec3* x, const PointData* pt)
+{
+ int index = umul24(blockIdx.x,BSIZE)+threadIdx.x;
+ if (index < size) {
+ if (!pt[index].bAccepted) {
+ x[index] = pt[index].oldpos;
+ }
+ }
+}
+
+
+
+//////////////////////
+// CPU-side methods //
+//////////////////////
+
+
+void Test2DAdapterCuda3f_computeTriangleNormal(unsigned int size, const void* x, void* tri)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_computeTriangleNormal_kernel<<< grid, threads >>>(size, (const CudaVec3*)x, (TriangleData*)tri);
+ mycudaDebugError("Test2DAdapterCuda3t_computeTriangleNormal_kernel");
+}
+
+void Test2DAdapterCuda3f_functionalGeom(unsigned int size, const void* x, void* tri)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_functionalGeom_kernel<<< grid, threads >>>(size, (const CudaVec3*)x, (TriangleData*)tri);
+ mycudaDebugError("Test2DAdapterCuda3t_functionalGeom_kernel");
+}
+
+void Test2DAdapterCuda3f_reduceStep(unsigned int size, void* x, void* pt, void* indices)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_reduceStep_kernel<<< grid, threads >>>(size, (CudaVec3*)x, (const PointData*)pt, (const Index*) indices);
+ mycudaDebugError("Test2DAdapterCuda3t_reduceStep_kernel");
+}
+
+void Test2DAdapterCuda3f_restoreUnchanged(unsigned int size, void* x, void* pt, void* indices)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_restoreUnchanged_kernel<<< grid, threads >>>(size, (CudaVec3*)x, (const PointData*)pt, (const Index*) indices);
+ mycudaDebugError("Test2DAdapterCuda3t_restoreUnchanged_kernel");
+}
+
+void Test2DAdapterCuda3f_smooth(unsigned int size, void* x, void* tri, void* pt, void* indices)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ //Test2DAdapterCuda3t_smoothLaplacian_kernel<<< grid, threads >>>(
+ // size, (CudaVec3*)x, (const TriangleData*)tri,
+ // (PointData*)pt, (const Index*) indices);
+ //mycudaDebugError("Test2DAdapterCuda3t_smoothLaplacian_kernel");
+ Test2DAdapterCuda3t_smoothOptimize_kernel<<< grid, threads >>>(
+ size, (CudaVec3*)x, (const TriangleData*)tri,
+ (PointData*)pt, (const Index*) indices);
+ mycudaDebugError("Test2DAdapterCuda3t_smoothOptimize_kernel");
+}
+
+void Test2DAdapterCuda3f_testAcceptable(unsigned int size, void* x, const void* tri, void* pt, void* indices, float tolerance)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_testAcceptable_kernel<<< grid, threads >>>(size, (CudaVec3*)x, (const TriangleData*)tri, (PointData*)pt, (const Index*) indices, tolerance);
+ mycudaDebugError("Test2DAdapterCuda3t_testAcceptable_kernel");
+}
+
+
+/// Specific to parallel version
+
+void Test2DAdapterCuda3f_prepareGradients(unsigned int size, const void* x,
+ void* tri)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_prepareGradients_kernel
+ <<< grid, threads >>>(
+ size, (CudaVec3*)x, (TriangleData*)tri);
+ mycudaDebugError("Test2DAdapterCuda3t_prepareGradients_kernel");
+}
+
+void Test2DAdapterCuda3f_smoothParallel(unsigned int size, void* x,
+ const void* tri, void* pt)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_smoothParallel_kernel
+ <<< grid, threads >>>(
+ size, (CudaVec3*)x, (const TriangleData*)tri,
+ (PointData*)pt);
+ mycudaDebugError("Test2DAdapterCuda3t_smoothParallel_kernel");
+}
+
+void Test2DAdapterCuda3f_reduceStepP(unsigned int size, void* x, void* pt)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_reduceStepP_kernel
+ <<< grid, threads >>>(
+ size, (CudaVec3*)x, (const PointData*)pt);
+ mycudaDebugError("Test2DAdapterCuda3t_reduceStepP_kernel");
+}
+
+void Test2DAdapterCuda3f_testAcceptableP(unsigned int size, void* x,
+ const void* tri, void* pt, float tolerance)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_testAcceptableP_kernel
+ <<< grid, threads >>>(
+ size, (CudaVec3*)x, (const TriangleData*)tri,
+ (PointData*)pt, tolerance);
+ mycudaDebugError("Test2DAdapterCuda3t_testAcceptableP_kernel");
+}
+
+void Test2DAdapterCuda3f_restoreUnchangedP(unsigned int size, void* x,
+ void* pt)
+{
+ dim3 threads(BSIZE,1);
+ dim3 grid((size+BSIZE-1)/BSIZE,1);
+ Test2DAdapterCuda3t_restoreUnchangedP_kernel
+ <<< grid, threads >>>(
+ size, (CudaVec3*)x, (const PointData*)pt);
+ mycudaDebugError("Test2DAdapterCuda3t_restoreUnchangedP_kernel");
+}
+
+
+#if defined(__cplusplus) && CUDA_VERSION < 2000
+} // namespace cuda
+} // namespace gpu
+} // namespace sofa
+#endif
diff --git a/src/Shell/controller/CudaTest2DAdapter.h b/src/Shell/controller/CudaTest2DAdapter.h
new file mode 100644
index 0000000..048c1e7
--- /dev/null
+++ b/src/Shell/controller/CudaTest2DAdapter.h
@@ -0,0 +1,110 @@
+#ifndef CUDA_TEST2DADAPTER_H
+#define CUDA_TEST2DADAPTER_H
+
+#include
+#include
+
+#include
+#include
+
+
+namespace sofa
+{
+
+namespace gpu
+{
+
+namespace cuda
+{
+
+template
+class CudaKernelsTest2DAdapter;
+
+} // namespace cuda
+
+} // namespace gpu
+
+namespace component
+{
+
+namespace controller
+{
+
+/**
+ * @brief Internal data for Test2DAdapter specialized for Cuda version.
+ */
+template<>
+class Test2DAdapterData< gpu::cuda::CudaVec3fTypes >
+{
+public:
+ typedef gpu::cuda::CudaVec3fTypes DataTypes;
+ typedef Test2DAdapter Main;
+ typedef DataTypes::Coord Coord;
+ typedef DataTypes::Deriv Deriv;
+ typedef Coord::value_type Real;
+
+ typedef sofa::component::topology::TriangleSetTopologyContainer::TriangleID Index;
+
+ // Graph colours (of independent sets)
+ type::vector< sofa::gpu::cuda::CudaVector > colours;
+
+ struct TriangleData {
+ type::fixed_array nodes;
+ Coord normal;
+ Real functional;
+ type::fixed_array gradient;
+ };
+
+ struct PointData {
+
+ bool bBoundary; /// Is the pont on a border?
+
+ unsigned int nNeighboursPt; /// Number of points in N1-ring.
+ const Index *neighboursPt; /// List of points in N1-ring.
+
+ unsigned int nNeighboursTri; /// Number of triangles in N1-ring.
+ const Index *neighboursTri; /// List of triangles in N1-ring.
+
+ bool bAccepted; /// New position has been accepted in current step.
+ Coord oldpos; /// Temporary holder for original point position.
+ Real oldworst;
+ Real newworst;
+
+ Index mintri; // Triangle from N1-ring with smallest functional
+ Deriv grad; // Gradient for mintri
+ };
+
+ struct PointDataHost {
+ sofa::gpu::cuda::CudaVector neighboursPt;
+ sofa::gpu::cuda::CudaVector neighboursTri;
+ };
+
+ sofa::gpu::cuda::CudaVector triangles;
+ sofa::gpu::cuda::CudaVector points;
+ sofa::gpu::cuda::CudaVector pointsHost;
+
+};
+
+template<>
+class Test2DAdapterData< gpu::cuda::CudaVec3fTypes >;
+
+template<>
+void Test2DAdapter< gpu::cuda::CudaVec3fTypes >::onEndAnimationStep(const double dt);
+
+template<>
+void Test2DAdapter< gpu::cuda::CudaVec3fTypes >::smoothLinear();
+
+template<>
+void Test2DAdapter< gpu::cuda::CudaVec3fTypes >::smoothParallel();
+
+template<>
+void Test2DAdapter< gpu::cuda::CudaVec3fTypes >::colourGraph();
+
+
+} // namespace controller
+
+} // namespace component
+
+} // namespace sofa
+
+#endif // CUDA_TEST2DADAPTER_H
diff --git a/src/Shell/controller/CudaTest2DAdapter.inl b/src/Shell/controller/CudaTest2DAdapter.inl
new file mode 100644
index 0000000..05223ae
--- /dev/null
+++ b/src/Shell/controller/CudaTest2DAdapter.inl
@@ -0,0 +1,390 @@
+#ifndef CUDA_TEST2DADAPTER_INL
+#define CUDA_TEST2DADAPTER_INL
+
+#include
+#include
+#include
+#include
+
+extern "C"
+{
+void Test2DAdapterCuda3f_computeTriangleNormal(unsigned int size, const void* x, const void* tri);
+void Test2DAdapterCuda3f_functionalGeom(unsigned int size, const void* x, void* tri);
+void Test2DAdapterCuda3f_reduceStep(unsigned int size, void* x, void* pt, void* indices);
+void Test2DAdapterCuda3f_restoreUnchanged(unsigned int size, void* x, void* pt, void* indices);
+void Test2DAdapterCuda3f_smooth(unsigned int size, void* x, void* tri, void* pt, void* indices);
+void Test2DAdapterCuda3f_testAcceptable(unsigned int size, void* x, const void* tri, void* pt, void* indices, float tolerance);
+///// parallel
+void Test2DAdapterCuda3f_prepareGradients(unsigned int size, const void* x, void* tri);
+void Test2DAdapterCuda3f_smoothParallel(unsigned int size, void* x, const void* tri, void* pt);
+void Test2DAdapterCuda3f_reduceStepP(unsigned int size, void* x, void* pt);
+void Test2DAdapterCuda3f_testAcceptableP(unsigned int size, void* x, const void* tri, void* pt, float tolerance);
+void Test2DAdapterCuda3f_restoreUnchangedP(unsigned int size, void* x, void* pt);
+//#ifdef SOFA_GPU_CUDA_DOUBLE
+//void Test2DAdapterCuda3d_computeTriangleNormal(const void* x, const void* n);
+//#endif
+}// extern "C"
+
+
+namespace sofa
+{
+
+namespace component
+{
+
+namespace controller
+{
+
+template<>
+void Test2DAdapter< gpu::cuda::CudaVec3fTypes >::onEndAnimationStep(const double /*dt*/)
+{
+ using namespace sofa::gpu::cuda;
+
+ std::cout << "GPU step\n";
+
+ if ((m_container == NULL) || (m_state == NULL))
+ return;
+
+ Index nTriangles = m_container->getNbTriangles();
+ if (nTriangles == 0)
+ return;
+
+ Data* datax = m_state->write(sofa::core::VecCoordId::position());
+ VecCoord& x = *datax->beginEdit();
+
+ //////////////////////////
+ // Initialization
+ if (data.colours.size() == 0) {
+ // Do the initialization
+ // TODO: move this elsewhere
+
+ // Colour the graph of vertices/edges to get the independent sets.
+ colourGraph();
+ // TODO: use a better container for the sets?
+
+ // Initialize triangle data
+ data.triangles.resize(nTriangles);
+ for (Index i=0; igetTriangle(i);
+ }
+
+ // Initialize point data
+ Index nPoints = x.size();
+ data.points.resize(nPoints);
+ data.pointsHost.resize(nPoints);
+ for (Index v=0; vgetEdgesAroundVertex(v);
+ data.pointsHost[v].neighboursPt.resize(N1e.size());
+
+ for (Index ie=0; iegetEdge(N1e[ie]);
+ for (int n=0; n<2; n++) {
+ if (e[n] != v) {
+ data.pointsHost[v].neighboursPt[ie] = e[n];
+ }
+ }
+ }
+
+ data.points[v].nNeighboursPt = N1e.size();
+ data.points[v].neighboursPt = reinterpret_cast
+ (data.pointsHost[v].neighboursPt.deviceRead());
+
+ // List of triangle neighbours
+ TrianglesAroundVertex N1 = m_container->getTrianglesAroundVertex(v);
+ data.pointsHost[v].neighboursTri.resize(N1.size());
+
+ for (Index it=0; it
+ (data.pointsHost[v].neighboursTri.deviceRead());
+ }
+ }
+ // End of Initialization
+ //////////////////////////
+
+ stepCounter++;
+
+ // Array sizes
+ //unsigned int szNormals = nTriangles * 3 * sizeof(Real);
+
+ // Compute initial metrics and normals
+ // TODO: kernel
+ //Real *normals_gpu;
+ //mycudaMalloc((void**)&normals_gpu, szNormals);
+
+ //Real *normals_buf;
+ //normals_buf = (Real*) malloc(szNormals);
+ //mycudaMemcpyDeviceToHost(normals_buf, normals_gpu, szNormals);
+
+ //vector normals(nTriangles);
+ //for (unsigned int i=0; i &functionals = *m_functionals.beginEdit();
+ functionals.resize(nTriangles);
+ for (Index i=0; iendEdit();
+}
+
+
+template<>
+void Test2DAdapter< gpu::cuda::CudaVec3fTypes >::smoothLinear()
+{
+ Index nTriangles = m_container->getNbTriangles();
+ if (nTriangles == 0)
+ return;
+
+ Data* datax = m_state->write(sofa::core::VecCoordId::position());
+ VecCoord& x = *datax->beginEdit();
+
+ // Persistent on gpu, we don't need to use this on host
+ void *triangles_gpu = data.triangles.deviceWrite();
+ void *points_gpu = data.points.deviceWrite();
+
+ // Compute initial normals
+ Test2DAdapterCuda3f_computeTriangleNormal(nTriangles, x.deviceRead(),
+ triangles_gpu);
+
+ // Compute initial values of the functional
+ Test2DAdapterCuda3f_functionalGeom(nTriangles, x.deviceRead(),
+ triangles_gpu);
+
+ //data.triangles.hostRead();
+ //std::cout << "m: ";
+ //for (unsigned int i=0; iendEdit();
+}
+
+template<>
+void Test2DAdapter< gpu::cuda::CudaVec3fTypes >::smoothParallel()
+{
+ Index nTriangles = m_container->getNbTriangles();
+ Index nPoints = m_container->getNbPoints();
+ if (nTriangles == 0)
+ return;
+
+ Data* datax = m_state->write(sofa::core::VecCoordId::position());
+ VecCoord& x = *datax->beginEdit();
+
+ // Persistent on gpu, we don't need to use this on host
+ void *triangles_gpu = data.triangles.deviceWrite();
+ void *points_gpu = data.points.deviceWrite();
+
+ // Compute initial normals
+ Test2DAdapterCuda3f_computeTriangleNormal(nTriangles, x.deviceRead(),
+ triangles_gpu);
+
+ // Compute initial values of the functional
+ Test2DAdapterCuda3f_functionalGeom(nTriangles, x.deviceRead(),
+ triangles_gpu);
+
+ // Prepare gradients
+ Test2DAdapterCuda3f_prepareGradients(nTriangles, x.deviceRead(),
+ triangles_gpu);
+
+ // Per vertex operations: gradient, assumed step
+ Test2DAdapterCuda3f_smoothParallel(nPoints, x.deviceWrite(),
+ triangles_gpu, points_gpu);
+
+ // Re-evaluate functional
+ Test2DAdapterCuda3f_functionalGeom(nTriangles, x.deviceRead(),
+ triangles_gpu);
+
+ // Test for acceptance
+ Test2DAdapterCuda3f_testAcceptableP(nPoints, x.deviceWrite(),
+ triangles_gpu, points_gpu, m_sigma.getValue());
+
+ data.points.hostRead();
+ for (unsigned int i=0; iendEdit();
+}
+
+template<>
+void Test2DAdapter< gpu::cuda::CudaVec3fTypes >::colourGraph()
+{
+ data.colours.clear();
+
+ int ncolours = 0;
+ type::vector c(m_container->getNbPoints(), -1);
+
+ for (Index v=0; (int)vgetNbPoints(); v++) {
+
+ if (pointInfo.getValue()[v].isBoundary() || pointInfo.getValue()[v].isFixed()) continue; // Skip boundary vertices
+
+ c[v] = 0;
+
+ EdgesAroundVertex N1e = m_container->getEdgesAroundVertex(v);
+ Index ie = 0;
+ while (ie < N1e.size()) {
+ Edge e = m_container->getEdge(N1e[ie]);
+ Index other = (e[0] == v) ? e[1] : e[0];
+ if (c[v] == c[other]) {
+ c[v]++;
+ ie = 0;
+ continue;
+ }
+ ie++;
+ }
+
+ if (c[v] >= ncolours) {
+ ncolours = c[v]+1;
+ data.colours.resize(ncolours);
+ }
+ data.colours[ c[v] ].push_back(v);
+ }
+
+ std::cout << "GC: " << data.colours.size() << " colours\n";
+}
+
+} // namespace controller
+
+} // namespace component
+
+} // namespace sofa
+
+#endif // SOFA_COMPONENT_CONTROLLER_TEST2DADAPTER_INL
diff --git a/src/Shell/controller/MeshChangedEvent.h b/src/Shell/controller/MeshChangedEvent.h
index ffdcf38..2a5aca5 100644
--- a/src/Shell/controller/MeshChangedEvent.h
+++ b/src/Shell/controller/MeshChangedEvent.h
@@ -32,7 +32,7 @@ namespace shell::objectmodel
using namespace sofa::type;
// Sent by MeshInterpolator when the mesh changes
-class MeshChangedEvent : public sofa::core::objectmodel::Event
+class SOFA_SHELL_API MeshChangedEvent : public sofa::core::objectmodel::Event
{
public:
diff --git a/src/Shell/controller/Test2DAdapter.cpp b/src/Shell/controller/Test2DAdapter.cpp
new file mode 100644
index 0000000..526c943
--- /dev/null
+++ b/src/Shell/controller/Test2DAdapter.cpp
@@ -0,0 +1,42 @@
+#include
+#include
+#include
+
+
+namespace sofa
+{
+
+namespace component
+{
+
+namespace controller
+{
+
+using namespace sofa::defaulttype;
+
+SOFA_DECL_CLASS(Test2DAdapter)
+
+// Register in the Factory
+int Test2DAdapterClass = core::RegisterObject("Adaptive mesh improvement component for 2D triangular meshes (for testing)")
+#ifdef SOFA_FLOAT
+.add< Test2DAdapter >(true) // default template
+#else
+.add< Test2DAdapter >(true) // default template
+# ifndef SOFA_DOUBLE
+.add< Test2DAdapter >()
+# endif
+#endif
+;
+
+#ifndef SOFA_FLOAT
+template class SOFA_SHELLS_API Test2DAdapter;
+#endif //SOFA_FLOAT
+#ifndef SOFA_DOUBLE
+template class SOFA_SHELLS_API Test2DAdapter;
+#endif //SOFA_DOUBLE
+
+} // namespace controller
+
+} // namespace component
+
+} // namespace sofa
diff --git a/src/Shell/controller/Test2DAdapter.h b/src/Shell/controller/Test2DAdapter.h
new file mode 100644
index 0000000..28f98bc
--- /dev/null
+++ b/src/Shell/controller/Test2DAdapter.h
@@ -0,0 +1,464 @@
+#ifndef SOFA_COMPONENT_CONTROLLER_TEST2DADAPTER_H
+#define SOFA_COMPONENT_CONTROLLER_TEST2DADAPTER_H
+
+#include
+
+#include
+#include
+
+#include