diff --git a/MbDCode/AtPointConstraintIJ.cpp b/MbDCode/AtPointConstraintIJ.cpp
index 9398f3b..b82209e 100644
--- a/MbDCode/AtPointConstraintIJ.cpp
+++ b/MbDCode/AtPointConstraintIJ.cpp
@@ -31,6 +31,8 @@ void AtPointConstraintIJ::initriIeJeO()
 
 void MbD::AtPointConstraintIJ::postInput()
 {
+	riIeJeO->postInput();
+	Constraint::postInput();
 }
 
 void MbD::AtPointConstraintIJ::calcPostDynCorrectorIteration()
@@ -43,3 +45,8 @@ void MbD::AtPointConstraintIJ::prePosIC()
 	riIeJeO->prePosIC();
 	Constraint::prePosIC();
 }
+
+MbD::ConstraintType MbD::AtPointConstraintIJ::type()
+{
+	return MbD::displacement;
+}
diff --git a/MbDCode/AtPointConstraintIJ.h b/MbDCode/AtPointConstraintIJ.h
index 1228a58..f5c52ce 100644
--- a/MbDCode/AtPointConstraintIJ.h
+++ b/MbDCode/AtPointConstraintIJ.h
@@ -17,6 +17,7 @@ namespace MbD {
         void postInput() override;
         void calcPostDynCorrectorIteration() override;
         void prePosIC() override;
+        MbD::ConstraintType type() override;
 
         int axis;
         std::shared_ptr<DispCompIecJecO> riIeJeO;
diff --git a/MbDCode/AtPointConstraintIqctJqc.cpp b/MbDCode/AtPointConstraintIqctJqc.cpp
index 6387e99..d12feea 100644
--- a/MbDCode/AtPointConstraintIqctJqc.cpp
+++ b/MbDCode/AtPointConstraintIqctJqc.cpp
@@ -27,3 +27,8 @@ void MbD::AtPointConstraintIqctJqc::calcPostDynCorrectorIteration()
 	ppGpEIpEI = std::static_pointer_cast<DispCompIeqctJeqcO>(riIeJeO)->ppriIeJeOpEIpEI;
 	AtPointConstraintIqcJqc::calcPostDynCorrectorIteration();
 }
+
+MbD::ConstraintType MbD::AtPointConstraintIqctJqc::type()
+{
+	return MbD::essential;
+}
diff --git a/MbDCode/AtPointConstraintIqctJqc.h b/MbDCode/AtPointConstraintIqctJqc.h
index b8432aa..f786cf0 100644
--- a/MbDCode/AtPointConstraintIqctJqc.h
+++ b/MbDCode/AtPointConstraintIqctJqc.h
@@ -11,6 +11,7 @@ namespace MbD {
         void initializeGlobally() override;
         void initriIeJeO() override;
         void calcPostDynCorrectorIteration() override;
+        MbD::ConstraintType type() override;
 
         double pGpt;
         FRowDsptr ppGpEIpt;
diff --git a/MbDCode/Constraint.cpp b/MbDCode/Constraint.cpp
index 790c8f5..0a9f877 100644
--- a/MbDCode/Constraint.cpp
+++ b/MbDCode/Constraint.cpp
@@ -1,5 +1,5 @@
 #include "Constraint.h"
-//#include "Item.h"
+#include "enum.h"
 
 using namespace MbD;
 
@@ -20,6 +20,8 @@ void Constraint::initialize()
 
 void MbD::Constraint::postInput()
 {
+	lam = 0.0;
+	Item::postInput();
 }
 
 void Constraint::setOwner(Item* x)
@@ -38,3 +40,27 @@ void MbD::Constraint::prePosIC()
 	iG = -1;
 	Item::prePosIC();
 }
+
+void MbD::Constraint::fillEssenConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints)
+{
+	if (this->type() == MbD::essential) {
+		essenConstraints->push_back(sptr);
+	}
+}
+void MbD::Constraint::fillDispConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> dispConstraints)
+{
+	if (this->type() == MbD::displacement) {
+		dispConstraints->push_back(sptr);
+	}
+}
+void MbD::Constraint::fillPerpenConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> perpenConstraints)
+{
+	if (this->type() == MbD::displacement) {
+		perpenConstraints->push_back(sptr);
+	}
+}
+
+MbD::ConstraintType MbD::Constraint::type()
+{
+	return MbD::essential;
+}
diff --git a/MbDCode/Constraint.h b/MbDCode/Constraint.h
index c6ed8af..54d2c5f 100644
--- a/MbDCode/Constraint.h
+++ b/MbDCode/Constraint.h
@@ -1,6 +1,7 @@
 #pragma once
 #include <memory>
 
+#include "enum.h"
 #include "Item.h"
 
 namespace MbD {
@@ -15,6 +16,10 @@ namespace MbD {
 		void setOwner(Item* x);
 		Item* getOwner();
 		void prePosIC() override;
+		virtual void fillEssenConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints);
+		virtual void fillDispConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> dispConstraints);
+		virtual void fillPerpenConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> perpenConstraints);
+		virtual MbD::ConstraintType type();
 
 		
 		int iG;
diff --git a/MbDCode/DirectionCosineConstraintIJ.cpp b/MbDCode/DirectionCosineConstraintIJ.cpp
index 9ea51c7..9105670 100644
--- a/MbDCode/DirectionCosineConstraintIJ.cpp
+++ b/MbDCode/DirectionCosineConstraintIJ.cpp
@@ -32,6 +32,8 @@ void DirectionCosineConstraintIJ::initaAijIeJe()
 
 void MbD::DirectionCosineConstraintIJ::postInput()
 {
+	aAijIeJe->postInput();
+	Constraint::postInput();
 }
 
 void MbD::DirectionCosineConstraintIJ::calcPostDynCorrectorIteration()
@@ -44,3 +46,8 @@ void MbD::DirectionCosineConstraintIJ::prePosIC()
 	aAijIeJe->prePosIC();
 	Constraint::prePosIC();
 }
+
+MbD::ConstraintType MbD::DirectionCosineConstraintIJ::type()
+{
+	return MbD::perpendicular;
+}
diff --git a/MbDCode/DirectionCosineConstraintIJ.h b/MbDCode/DirectionCosineConstraintIJ.h
index 12470ec..18e9ab1 100644
--- a/MbDCode/DirectionCosineConstraintIJ.h
+++ b/MbDCode/DirectionCosineConstraintIJ.h
@@ -17,6 +17,7 @@ namespace MbD {
         void postInput() override;
         void calcPostDynCorrectorIteration() override;
         void prePosIC() override;
+        MbD::ConstraintType type() override;
 
         int axisI, axisJ;
         std::shared_ptr<DirectionCosineIecJec> aAijIeJe;
diff --git a/MbDCode/DirectionCosineConstraintIqctJqc.cpp b/MbDCode/DirectionCosineConstraintIqctJqc.cpp
index a15cdbf..b27bc59 100644
--- a/MbDCode/DirectionCosineConstraintIqctJqc.cpp
+++ b/MbDCode/DirectionCosineConstraintIqctJqc.cpp
@@ -13,3 +13,8 @@ void DirectionCosineConstraintIqctJqc::initaAijIeJe()
 {
 	aAijIeJe = CREATE<DirectionCosineIeqctJeqc>::With(frmI, frmJ, axisI, axisJ);
 }
+
+MbD::ConstraintType MbD::DirectionCosineConstraintIqctJqc::type()
+{
+	return MbD::essential;
+}
diff --git a/MbDCode/DirectionCosineConstraintIqctJqc.h b/MbDCode/DirectionCosineConstraintIqctJqc.h
index 9ed5c7b..38e2fa1 100644
--- a/MbDCode/DirectionCosineConstraintIqctJqc.h
+++ b/MbDCode/DirectionCosineConstraintIqctJqc.h
@@ -9,6 +9,7 @@ namespace MbD {
 	public:
 		DirectionCosineConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
 		void initaAijIeJe() override;
+		MbD::ConstraintType type() override;
 
 		double pGpt;
 		FRowDsptr ppGpEIpt;
diff --git a/MbDCode/DispCompIeqctJeqcO.cpp b/MbDCode/DispCompIeqctJeqcO.cpp
index 024ed44..5f96018 100644
--- a/MbDCode/DispCompIeqctJeqcO.cpp
+++ b/MbDCode/DispCompIeqctJeqcO.cpp
@@ -12,12 +12,12 @@ MbD::DispCompIeqctJeqcO::DispCompIeqctJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, in
 void MbD::DispCompIeqctJeqcO::initializeGlobally()
 {
 	//ToDo: Check why not using super classes.
-	ppriIeJeOpEJpEJ = std::static_pointer_cast<EndFrameqc>(frmJ)->ppriOeOpEpE(axis);
+	ppriIeJeOpEJpEJ = std::static_pointer_cast<EndFrameqct>(frmJ)->ppriOeOpEpE(axis);
 }
 
 void MbD::DispCompIeqctJeqcO::calcPostDynCorrectorIteration()
 {
 	//"ppriIeJeOpEIpEI is not a constant now."
 	DispCompIeqcJeqcO::calcPostDynCorrectorIteration();
-	ppriIeJeOpEIpEI = std::static_pointer_cast<EndFrameqc>(frmI)->ppriOeOpEpE(axis)->negated();
+	ppriIeJeOpEIpEI = std::static_pointer_cast<EndFrameqct>(frmI)->ppriOeOpEpE(axis)->negated();
 }
diff --git a/MbDCode/EndFrameqc.cpp b/MbDCode/EndFrameqc.cpp
index b1cdc48..f1c607d 100644
--- a/MbDCode/EndFrameqc.cpp
+++ b/MbDCode/EndFrameqc.cpp
@@ -35,6 +35,7 @@ void EndFrameqc::initEndFrameqct()
 	endFrameqct->pprOeOpEpE = pprOeOpEpE;
 	endFrameqct->pAOepE = pAOepE;
 	endFrameqct->ppAOepEpE = ppAOepEpE;
+	endFrameqct->setMarkerFrame(markerFrame);
 }
 
 FMatFColDsptr MbD::EndFrameqc::ppAjOepEpE(int jj)
diff --git a/MbDCode/EndFrameqct.cpp b/MbDCode/EndFrameqct.cpp
index 3ae7cd5..5ae3a7b 100644
--- a/MbDCode/EndFrameqct.cpp
+++ b/MbDCode/EndFrameqct.cpp
@@ -92,6 +92,9 @@ void MbD::EndFrameqct::initppPhiThePsiptptBlks()
 
 void MbD::EndFrameqct::postInput()
 {
+	this->evalrmem();
+	this->evalAme();
+	Item::postInput();
 }
 
 void MbD::EndFrameqct::calcPostDynCorrectorIteration()
diff --git a/MbDCode/EulerConstraint.cpp b/MbDCode/EulerConstraint.cpp
index 132e21f..10a5947 100644
--- a/MbDCode/EulerConstraint.cpp
+++ b/MbDCode/EulerConstraint.cpp
@@ -20,7 +20,7 @@ void EulerConstraint::initialize()
 
 void MbD::EulerConstraint::calcPostDynCorrectorIteration()
 {
-	auto qE = static_cast<PartFrame*>(owner)->qE;
+	auto& qE = static_cast<PartFrame*>(owner)->qE;
 	aG = qE->sumOfSquares() - 1.0;
 	for (size_t i = 0; i < 4; i++)
 	{
diff --git a/MbDCode/ForceTorqueItem.cpp b/MbDCode/ForceTorqueItem.cpp
new file mode 100644
index 0000000..9049fa9
--- /dev/null
+++ b/MbDCode/ForceTorqueItem.cpp
@@ -0,0 +1 @@
+#include "ForceTorqueItem.h"
diff --git a/MbDCode/ForceTorqueItem.h b/MbDCode/ForceTorqueItem.h
new file mode 100644
index 0000000..3518953
--- /dev/null
+++ b/MbDCode/ForceTorqueItem.h
@@ -0,0 +1,13 @@
+#pragma once
+
+#include "Item.h"
+
+namespace MbD {
+    class ForceTorqueItem : public Item
+    {
+        //
+    public:
+
+    };
+}
+
diff --git a/MbDCode/FullRow.h b/MbDCode/FullRow.h
index 8bb0cb2..a55f857 100644
--- a/MbDCode/FullRow.h
+++ b/MbDCode/FullRow.h
@@ -72,7 +72,13 @@ namespace MbD {
 	template<typename T>
 	inline std::shared_ptr<FullRow<T>> FullRow<T>::timesTransposeFullMatrix(std::shared_ptr<FullMatrix<T>> fullMat)
 	{
-		return std::shared_ptr<FullRow<T>>();
+		//"a*bT = a(1,j)b(k,j)"
+		size_t ncol = fullMat->nRow();
+		auto answer = std::make_shared<FullRow<T>>(ncol);
+		for (size_t k = 0; k < ncol; k++) {
+			answer->at(k) = this->dot(fullMat->at(k));
+		}
+		return answer;
 	}
 	template<typename T>
 	inline void FullRow<T>::equalSelfPlusFullRowTimes(std::shared_ptr<FullRow<T>> fullRow, double factor)
diff --git a/MbDCode/Item.cpp b/MbDCode/Item.cpp
index 885d821..e3f87da 100644
--- a/MbDCode/Item.cpp
+++ b/MbDCode/Item.cpp
@@ -57,6 +57,10 @@ void MbD::Item::fillPosKineJacob(FMatDsptr mat)
 {
 }
 
+void MbD::Item::fillEssenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints)
+{
+}
+
 void MbD::Item::constraintsReport()
 {
 }
diff --git a/MbDCode/Item.h b/MbDCode/Item.h
index 1e48da2..1ca3f51 100644
--- a/MbDCode/Item.h
+++ b/MbDCode/Item.h
@@ -7,6 +7,8 @@
 
 namespace MbD {
 
+	class Constraint;
+
 	class Item
 	{
 		//name
@@ -31,6 +33,7 @@ namespace MbD {
 		virtual void removeRedundantConstraints(std::shared_ptr<std::vector<int>> redundantEqnNos);
 		virtual void fillPosKineError(FColDsptr col);
 		virtual void fillPosKineJacob(FMatDsptr mat);
+		virtual void fillEssenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints);
 
 		void setName(std::string& str);
 		const std::string& getName() const;
diff --git a/MbDCode/Joint.cpp b/MbDCode/Joint.cpp
index 5511abd..df85836 100644
--- a/MbDCode/Joint.cpp
+++ b/MbDCode/Joint.cpp
@@ -34,16 +34,23 @@ void Joint::initializeLocally()
 			frmI = frmIqc->endFrameqct;
 		}
 	}
-	std::for_each(constraints->begin(), constraints->end(), [](const auto& constraint) { constraint->initializeLocally(); });
+	constraintsDo([](const auto& constraint) { constraint->initializeLocally(); });
 }
 
 void Joint::initializeGlobally()
 {
-	std::for_each(constraints->begin(), constraints->end(), [](const auto& constraint) { constraint->initializeGlobally(); });
+	constraintsDo([](const auto& constraint) { constraint->initializeGlobally(); });
+}
+
+void MbD::Joint::constraintsDo(const std::function<void(std::shared_ptr<Constraint>)>& f)
+{
+	std::for_each(constraints->begin(), constraints->end(), f);
 }
 
 void MbD::Joint::postInput()
 {
+	constraintsDo([](const auto& constraint) { constraint->postInput(); });
+
 }
 
 void MbD::Joint::addConstraint(std::shared_ptr<Constraint> con)
@@ -54,5 +61,20 @@ void MbD::Joint::addConstraint(std::shared_ptr<Constraint> con)
 
 void MbD::Joint::prePosIC()
 {
-	std::for_each(constraints->begin(), constraints->end(), [](const auto& constraint) { constraint->prePosIC(); });
+	constraintsDo([](const auto& constraint) { constraint->prePosIC(); });
+}
+
+void MbD::Joint::fillEssenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints)
+{
+	constraintsDo([&](const auto& con) { con->fillEssenConstraints(con, essenConstraints); });
+}
+
+void MbD::Joint::fillDispConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> dispConstraints)
+{
+	constraintsDo([&](const auto& con) { con->fillDispConstraints(con, dispConstraints); });
+}
+
+void MbD::Joint::fillPerpenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> perpenConstraints)
+{
+	constraintsDo([&](const auto& con) { con->fillPerpenConstraints(con, perpenConstraints); });
 }
diff --git a/MbDCode/Joint.h b/MbDCode/Joint.h
index 7143c66..0b8261c 100644
--- a/MbDCode/Joint.h
+++ b/MbDCode/Joint.h
@@ -1,6 +1,7 @@
 #pragma once
 #include <memory>
 #include <vector>
+#include <functional>
 
 #include "Item.h"
 #include "EndFramec.h"
@@ -18,9 +19,13 @@ namespace MbD {
         virtual void connectsItoJ(EndFrmcptr frmI, EndFrmcptr frmJ);
         void initializeLocally() override;
         void initializeGlobally() override;
+        void constraintsDo(const std::function <void(std::shared_ptr<Constraint>)>& f);
         void postInput() override;
         void addConstraint(std::shared_ptr<Constraint> con);
         void prePosIC() override;
+        void fillEssenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints) override;
+        virtual void fillDispConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> dispConstraints);
+        virtual void fillPerpenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> perpenConstraints);
 
         EndFrmcptr frmI;
         EndFrmcptr frmJ;
diff --git a/MbDCode/MarkerFrame.cpp b/MbDCode/MarkerFrame.cpp
index dfcfda6..00492c8 100644
--- a/MbDCode/MarkerFrame.cpp
+++ b/MbDCode/MarkerFrame.cpp
@@ -48,6 +48,8 @@ void MarkerFrame::initializeGlobally()
 
 void MbD::MarkerFrame::postInput()
 {
+	Item::postInput();
+	std::for_each(endFrames->begin(), endFrames->end(), [](const auto& endFrame) { endFrame->postInput(); });
 }
 
 void MbD::MarkerFrame::calcPostDynCorrectorIteration()
diff --git a/MbDCode/MbDCode.aps b/MbDCode/MbDCode.aps
new file mode 100644
index 0000000..39204fa
Binary files /dev/null and b/MbDCode/MbDCode.aps differ
diff --git a/MbDCode/MbDCode.rc b/MbDCode/MbDCode.rc
new file mode 100644
index 0000000..44acb21
Binary files /dev/null and b/MbDCode/MbDCode.rc differ
diff --git a/MbDCode/MbDCode.vcxproj b/MbDCode/MbDCode.vcxproj
index 4623cae..4ed474d 100644
--- a/MbDCode/MbDCode.vcxproj
+++ b/MbDCode/MbDCode.vcxproj
@@ -171,6 +171,7 @@
     <ClCompile Include="EulerParameters.cpp" />
     <ClCompile Include="EulerParametersDot.cpp" />
     <ClCompile Include="ExpressionX.cpp" />
+    <ClCompile Include="ForceTorqueItem.cpp" />
     <ClCompile Include="FullColumn.cpp" />
     <ClCompile Include="FullMatrix.cpp" />
     <ClCompile Include="FullRow.cpp" />
@@ -208,6 +209,7 @@
     <ClCompile Include="PrescribedMotion.cpp" />
     <ClCompile Include="Product.cpp" />
     <ClCompile Include="QuasiIntegrator.cpp" />
+    <ClCompile Include="RedundantConstraint.cpp" />
     <ClCompile Include="RevoluteJoint.cpp" />
     <ClCompile Include="RowTypeMatrix.cpp" />
     <ClCompile Include="ScalarNewtonRaphson.cpp" />
@@ -272,12 +274,14 @@
     <ClInclude Include="EndFramec.h" />
     <ClInclude Include="EndFrameqc.h" />
     <ClInclude Include="EndFrameqct.h" />
+    <ClInclude Include="enum.h" />
     <ClInclude Include="EulerAngleszxz.h" />
     <ClInclude Include="EulerArray.h" />
     <ClInclude Include="EulerConstraint.h" />
     <ClInclude Include="EulerParameters.h" />
     <ClInclude Include="EulerParametersDot.h" />
     <ClInclude Include="ExpressionX.h" />
+    <ClInclude Include="ForceTorqueItem.h" />
     <ClInclude Include="FullColumn.h" />
     <ClInclude Include="FullMatrix.h" />
     <ClInclude Include="FullRow.h" />
@@ -315,6 +319,8 @@
     <ClInclude Include="PrescribedMotion.h" />
     <ClInclude Include="Product.h" />
     <ClInclude Include="QuasiIntegrator.h" />
+    <ClInclude Include="RedundantConstraint.h" />
+    <ClInclude Include="resource.h" />
     <ClInclude Include="RevoluteJoint.h" />
     <ClInclude Include="RowTypeMatrix.h" />
     <ClInclude Include="ScalarNewtonRaphson.h" />
@@ -341,6 +347,9 @@
     <ClInclude Include="VectorNewtonRaphson.h" />
     <ClInclude Include="ZRotation.h" />
   </ItemGroup>
+  <ItemGroup>
+    <ResourceCompile Include="MbDCode.rc" />
+  </ItemGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">
   </ImportGroup>
diff --git a/MbDCode/MbDCode.vcxproj.filters b/MbDCode/MbDCode.vcxproj.filters
index c9cd1a1..950e479 100644
--- a/MbDCode/MbDCode.vcxproj.filters
+++ b/MbDCode/MbDCode.vcxproj.filters
@@ -330,6 +330,12 @@
     <ClCompile Include="GEFullMatParPv.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="ForceTorqueItem.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="RedundantConstraint.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="Array.h">
@@ -650,5 +656,22 @@
     <ClInclude Include="GEFullMatParPv.h">
       <Filter>Header Files</Filter>
     </ClInclude>
+    <ClInclude Include="ForceTorqueItem.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="resource.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="enum.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="RedundantConstraint.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+  </ItemGroup>
+  <ItemGroup>
+    <ResourceCompile Include="MbDCode.rc">
+      <Filter>Resource Files</Filter>
+    </ResourceCompile>
   </ItemGroup>
 </Project>
\ No newline at end of file
diff --git a/MbDCode/NewtonRaphson.cpp b/MbDCode/NewtonRaphson.cpp
index c6a96a2..c742e11 100644
--- a/MbDCode/NewtonRaphson.cpp
+++ b/MbDCode/NewtonRaphson.cpp
@@ -7,9 +7,9 @@ using namespace MbD;
 
 void NewtonRaphson::initialize()
 {
-	dxNorms = std::shared_ptr<std::vector<double>>();
+	dxNorms = std::make_shared<std::vector<double>>();
 	dxTol = 4 * std::numeric_limits<double>::epsilon();
-	yNorms = std::shared_ptr<std::vector<double>>();
+	yNorms = std::make_shared<std::vector<double>>();
 	yNormTol = 1.0e-30;
 	iterMax = 100;
 	twoAlp = 2.0e-4;
diff --git a/MbDCode/Part.cpp b/MbDCode/Part.cpp
index 0723082..24969be 100644
--- a/MbDCode/Part.cpp
+++ b/MbDCode/Part.cpp
@@ -76,6 +76,8 @@ void Part::asFixed()
 
 void MbD::Part::postInput()
 {
+	partFrame->postInput();
+	Item::postInput();
 }
 
 void MbD::Part::calcPostDynCorrectorIteration()
@@ -97,3 +99,8 @@ void MbD::Part::iqE(int eqnNo)
 	partFrame->iqE = eqnNo;
 
 }
+
+void MbD::Part::fillEssenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints)
+{
+	partFrame->fillEssenConstraints(essenConstraints);
+}
diff --git a/MbDCode/Part.h b/MbDCode/Part.h
index 5620e49..0c2134d 100644
--- a/MbDCode/Part.h
+++ b/MbDCode/Part.h
@@ -36,6 +36,7 @@ namespace MbD {
 		void prePosIC() override;
 		void iqX(int eqnNo);
 		void iqE(int eqnNo);
+		void fillEssenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints) override;
 
 		int ipX = -1; 
 		int ipE = -1; 
diff --git a/MbDCode/PartFrame.cpp b/MbDCode/PartFrame.cpp
index 365039b..6d6ad93 100644
--- a/MbDCode/PartFrame.cpp
+++ b/MbDCode/PartFrame.cpp
@@ -29,14 +29,14 @@ void PartFrame::initializeLocally()
 {
 	std::for_each(markerFrames->begin(), markerFrames->end(), [](const auto& markerFrame) { markerFrame->initializeLocally(); });
 	aGeu->initializeLocally();
-	std::for_each(aGabs->begin(), aGabs->end(), [](const auto& aGab) { aGab->initializeLocally(); });
+	aGabsDo([](const auto& aGab) { aGab->initializeLocally(); });
 }
 
 void PartFrame::initializeGlobally()
 {
 	std::for_each(markerFrames->begin(), markerFrames->end(), [](const auto& markerFrame) { markerFrame->initializeGlobally(); });
 	aGeu->initializeGlobally();
-	std::for_each(aGabs->begin(), aGabs->end(), [](const auto& aGab) { aGab->initializeGlobally(); });
+	aGabsDo([](const auto& aGab) { aGab->initializeGlobally(); });
 }
 void PartFrame::setqX(FColDsptr x) {
 	qX->copyFrom(x);
@@ -90,6 +90,11 @@ EndFrmcptr PartFrame::endFrame(std::string name)
 	return (*match)->endFrames->at(0);
 }
 
+void MbD::PartFrame::aGabsDo(const std::function<void(std::shared_ptr<AbsConstraint>)>& f)
+{
+	std::for_each(aGabs->begin(), aGabs->end(), f);
+}
+
 void MbD::PartFrame::prePosIC()
 {
 //iqX = -1;
@@ -115,6 +120,12 @@ FColFMatDsptr MbD::PartFrame::pAOppE()
 	return qE->pApE;
 }
 
+void MbD::PartFrame::fillEssenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints)
+{
+	aGeu->fillEssenConstraints(aGeu, essenConstraints);
+	aGabsDo([&](const auto& con) { con->fillEssenConstraints(con, essenConstraints); });
+}
+
 void PartFrame::asFixed()
 {
 	for (size_t i = 0; i < 6; i++) {
@@ -126,6 +137,12 @@ void PartFrame::asFixed()
 
 void MbD::PartFrame::postInput()
 {
+	//qXddot = std::make_shared<FullColumn<double>>(3, 0.0);
+	//qEddot = std::make_shared<FullColumn<double>>(4, 0.0);
+	Item::postInput();
+	std::for_each(markerFrames->begin(), markerFrames->end(), [](const auto& markerFrame) { markerFrame->postInput(); });
+	aGeu->postInput();
+	aGabsDo([](const auto& aGab) { aGab->postInput(); });
 }
 
 void MbD::PartFrame::calcPostDynCorrectorIteration()
diff --git a/MbDCode/PartFrame.h b/MbDCode/PartFrame.h
index b751e39..423e1bc 100644
--- a/MbDCode/PartFrame.h
+++ b/MbDCode/PartFrame.h
@@ -1,6 +1,7 @@
 #pragma once
 #include <memory>
 #include <vector>
+#include <functional>
 
 #include "CartesianFrame.h"
 #include "EndFramec.h"
@@ -42,11 +43,13 @@ namespace MbD {
 		Part* getPart();
 		void addMarkerFrame(std::shared_ptr<MarkerFrame> x);
 		EndFrmcptr endFrame(std::string name);
+		void aGabsDo(const std::function <void(std::shared_ptr<AbsConstraint>)>& f);
 
 		void prePosIC() override;
 		FColDsptr rOpO();
 		FMatDsptr aAOp();
 		FColFMatDsptr pAOppE();
+		void fillEssenConstraints(std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essenConstraints) override;
 
 		Part* part = nullptr;
 		int iqX = -1;
diff --git a/MbDCode/PosICNewtonRaphson.cpp b/MbDCode/PosICNewtonRaphson.cpp
index 73854c7..ee9d164 100644
--- a/MbDCode/PosICNewtonRaphson.cpp
+++ b/MbDCode/PosICNewtonRaphson.cpp
@@ -70,7 +70,7 @@ void MbD::PosICNewtonRaphson::assignEquationNumbers()
 		con->iG = eqnNo;
 		eqnNo = eqnNo + 1;
 	}
-	auto n = eqnNo;
+	n = eqnNo;
 	auto limits = { lastEssenConEqnNo, lastDispConEqnNo, n };
 	pivotRowLimits = std::make_shared<Array<int>>(limits);
 }
diff --git a/MbDCode/RedundantConstraint.cpp b/MbDCode/RedundantConstraint.cpp
new file mode 100644
index 0000000..49460e2
--- /dev/null
+++ b/MbDCode/RedundantConstraint.cpp
@@ -0,0 +1 @@
+#include "RedundantConstraint.h"
diff --git a/MbDCode/RedundantConstraint.h b/MbDCode/RedundantConstraint.h
new file mode 100644
index 0000000..9c25ddf
--- /dev/null
+++ b/MbDCode/RedundantConstraint.h
@@ -0,0 +1,13 @@
+#pragma once
+
+#include "Item.h"
+
+namespace MbD {
+    class RedundantConstraint : public Item
+    {
+        //
+    public:
+        std::shared_ptr<Constraint> constraint;
+    };
+}
+
diff --git a/MbDCode/System.cpp b/MbDCode/System.cpp
index b9157b9..82648bd 100644
--- a/MbDCode/System.cpp
+++ b/MbDCode/System.cpp
@@ -3,6 +3,7 @@
 #include "System.h"
 #include "Part.h"
 #include "Joint.h"
+#include "ForceTorqueItem.h"
 #include "SystemSolver.h"
 #include "Time.h"
 #include "CREATE.h"
@@ -22,6 +23,7 @@ void MbD::System::initialize()
 	time = CREATE<Time>::With();
 	parts = std::make_shared<std::vector<std::shared_ptr<Part>>>();
 	jointsMotions = std::make_shared<std::vector<std::shared_ptr<Joint>>>();
+	forcesTorques = std::make_shared<std::vector<std::shared_ptr<ForceTorqueItem>>>();
 	systemSolver = std::make_shared<SystemSolver>(this);
 }
 
@@ -39,7 +41,7 @@ void System::runKINEMATICS()
 		initializeGlobally();
 		if (!hasChanged) break;
 	}
-	postInput();
+	partsJointsMotionsForcesTorquesDo([&](std::shared_ptr<Item> item) { item->postInput(); });
 	outputInput();
 	systemSolver->runAllIC();
 	outputInitialConditions();
@@ -63,15 +65,13 @@ void System::initializeLocally()
 {
 	hasChanged = false;
 	time->value = systemSolver->tstart;
-	std::for_each(parts->begin(), parts->end(), [](const auto& part) { part->initializeLocally(); });
-	std::for_each(jointsMotions->begin(), jointsMotions->end(), [](const auto& joint) { joint->initializeLocally();	});
+	partsJointsMotionsDo([](const auto& item) { item->initializeLocally(); });
 	systemSolver->initializeLocally();
 }
 
 void System::initializeGlobally()
 {
-	std::for_each(parts->begin(), parts->end(), [](const auto& part) { part->initializeGlobally(); });
-	std::for_each(jointsMotions->begin(), jointsMotions->end(), [](const auto& joint) { joint->initializeGlobally();	});
+	partsJointsMotionsDo([](const auto& item) { item->initializeGlobally(); });
 	systemSolver->initializeGlobally();
 }
 
@@ -80,12 +80,24 @@ std::shared_ptr<std::vector<std::string>> MbD::System::discontinuitiesAtIC()
 	return std::shared_ptr<std::vector<std::string>>();
 }
 
+void MbD::System::jointsMotionsDo(const std::function<void(std::shared_ptr<Joint>)>& f)
+{
+	std::for_each(jointsMotions->begin(), jointsMotions->end(), f);
+}
+
 void MbD::System::partsJointsMotionsDo(const std::function<void(std::shared_ptr<Item>)>& f)
 {
 	std::for_each(parts->begin(), parts->end(), f);
 	std::for_each(jointsMotions->begin(), jointsMotions->end(), f);
 }
 
+void MbD::System::partsJointsMotionsForcesTorquesDo(const std::function<void(std::shared_ptr<Item>)>& f)
+{
+	std::for_each(parts->begin(), parts->end(), f);
+	std::for_each(jointsMotions->begin(), jointsMotions->end(), f);
+	std::for_each(forcesTorques->begin(), forcesTorques->end(), f);
+}
+
 void MbD::System::logString(std::string& str)
 {
 	std::cout << str << std::endl;
@@ -95,3 +107,24 @@ double MbD::System::mbdTimeValue()
 {
 	return time->getValue();
 }
+
+std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> MbD::System::essentialConstraints2()
+{
+	auto essenConstraints = std::make_shared<std::vector<std::shared_ptr<Constraint>>>();
+	this->partsJointsMotionsDo([&](const auto& item) { item->fillEssenConstraints(essenConstraints); });
+	return essenConstraints;
+}
+
+std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> MbD::System::displacementConstraints()
+{
+	auto dispConstraints = std::make_shared<std::vector<std::shared_ptr<Constraint>>>();
+	this->jointsMotionsDo([&](const auto& joint) { joint->fillDispConstraints(dispConstraints); });
+	return dispConstraints;
+}
+
+std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> MbD::System::perpendicularConstraints2()
+{
+	auto perpenConstraints = std::make_shared<std::vector<std::shared_ptr<Constraint>>>();
+	this->jointsMotionsDo([&](const auto& joint) { joint->fillPerpenConstraints(perpenConstraints); });
+	return perpenConstraints;
+}
diff --git a/MbDCode/System.h b/MbDCode/System.h
index 11d759e..53cc536 100644
--- a/MbDCode/System.h
+++ b/MbDCode/System.h
@@ -19,6 +19,8 @@ namespace MbD {
 	class Joint;
 	class SystemSolver;
 	class Time;
+	class Constraint;
+	class ForceTorqueItem;
 
 	class System : public Item
 	{
@@ -38,13 +40,19 @@ namespace MbD {
 		void outputInitialConditions();
 		void outputTimeSeries();
 		std::shared_ptr<std::vector<std::string>> discontinuitiesAtIC();
+		void jointsMotionsDo(const std::function <void(std::shared_ptr<Joint>)>& f);
 		void partsJointsMotionsDo(const std::function <void(std::shared_ptr<Item>)>& f);
+		void partsJointsMotionsForcesTorquesDo(const std::function <void(std::shared_ptr<Item>)>& f);
 		void logString(std::string& str);
 		double mbdTimeValue();
+		std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> essentialConstraints2();
+		std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> displacementConstraints();
+		std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> perpendicularConstraints2();
 
 
 		std::shared_ptr<std::vector<std::shared_ptr<Part>>> parts;
 		std::shared_ptr<std::vector<std::shared_ptr<Joint>>> jointsMotions;
+		std::shared_ptr<std::vector<std::shared_ptr<ForceTorqueItem>>> forcesTorques;
 		bool hasChanged = false;
 		std::shared_ptr<SystemSolver> systemSolver;
 		void addPart(std::shared_ptr<Part> part);
diff --git a/MbDCode/SystemSolver.cpp b/MbDCode/SystemSolver.cpp
index 0b5158e..6bdb392 100644
--- a/MbDCode/SystemSolver.cpp
+++ b/MbDCode/SystemSolver.cpp
@@ -97,20 +97,20 @@ void MbD::SystemSolver::logString(std::string& str)
 
 std::shared_ptr<std::vector<std::shared_ptr<Part>>> MbD::SystemSolver::parts()
 {
-	return std::shared_ptr<std::vector<std::shared_ptr<Part>>>();
+	return system->parts;
 }
 
 std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> MbD::SystemSolver::essentialConstraints2()
 {
-	return std::shared_ptr<std::vector<std::shared_ptr<Constraint>>>();
+	return system->essentialConstraints2();
 }
 
 std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> MbD::SystemSolver::displacementConstraints()
 {
-	return std::shared_ptr<std::vector<std::shared_ptr<Constraint>>>();
+	return system->displacementConstraints();
 }
 
 std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> MbD::SystemSolver::perpendicularConstraints2()
 {
-	return std::shared_ptr<std::vector<std::shared_ptr<Constraint>>>();
+	return system->perpendicularConstraints2();
 }
diff --git a/MbDCode/TranslationConstraintIJ.cpp b/MbDCode/TranslationConstraintIJ.cpp
index 303db62..a806bbe 100644
--- a/MbDCode/TranslationConstraintIJ.cpp
+++ b/MbDCode/TranslationConstraintIJ.cpp
@@ -31,6 +31,8 @@ void MbD::TranslationConstraintIJ::initriIeJeIe()
 
 void MbD::TranslationConstraintIJ::postInput()
 {
+    riIeJeIe->postInput();
+    Constraint::postInput();
 }
 
 void MbD::TranslationConstraintIJ::calcPostDynCorrectorIteration()
@@ -43,3 +45,8 @@ void MbD::TranslationConstraintIJ::prePosIC()
     riIeJeIe->prePosIC();
     Constraint::prePosIC();
 }
+
+MbD::ConstraintType MbD::TranslationConstraintIJ::type()
+{
+    return MbD::displacement;
+}
diff --git a/MbDCode/TranslationConstraintIJ.h b/MbDCode/TranslationConstraintIJ.h
index b0f9528..7960de9 100644
--- a/MbDCode/TranslationConstraintIJ.h
+++ b/MbDCode/TranslationConstraintIJ.h
@@ -16,6 +16,7 @@ namespace MbD {
         void postInput() override;
         void calcPostDynCorrectorIteration() override;
         void prePosIC()override;
+        MbD::ConstraintType type() override;
 
         int axisI;
         std::shared_ptr<DispCompIecJecKec> riIeJeIe;
diff --git a/MbDCode/TranslationConstraintIqctJqc.cpp b/MbDCode/TranslationConstraintIqctJqc.cpp
index dd7a9e1..58f0f43 100644
--- a/MbDCode/TranslationConstraintIqctJqc.cpp
+++ b/MbDCode/TranslationConstraintIqctJqc.cpp
@@ -13,3 +13,8 @@ void MbD::TranslationConstraintIqctJqc::initriIeJeIe()
 {
 	riIeJeIe = CREATE<DispCompIeqctJeqcKeqct>::With(frmI, frmJ, frmI, axisI);
 }
+
+MbD::ConstraintType MbD::TranslationConstraintIqctJqc::type()
+{
+	return MbD::essential;
+}
diff --git a/MbDCode/TranslationConstraintIqctJqc.h b/MbDCode/TranslationConstraintIqctJqc.h
index 8385d30..5e52f79 100644
--- a/MbDCode/TranslationConstraintIqctJqc.h
+++ b/MbDCode/TranslationConstraintIqctJqc.h
@@ -9,6 +9,7 @@ namespace MbD {
     public:
         TranslationConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi);
         void initriIeJeIe() override;
+        MbD::ConstraintType type() override;
 
         double pGpt;
         FRowDsptr ppGpXIpt;
diff --git a/MbDCode/enum.h b/MbDCode/enum.h
new file mode 100644
index 0000000..e867d9e
--- /dev/null
+++ b/MbDCode/enum.h
@@ -0,0 +1,4 @@
+#pragma once
+namespace MbD {
+	enum ConstraintType {essential, displacement, perpendicular};
+}
\ No newline at end of file
diff --git a/MbDCode/resource.h b/MbDCode/resource.h
new file mode 100644
index 0000000..0d5b956
--- /dev/null
+++ b/MbDCode/resource.h
@@ -0,0 +1,14 @@
+//{{NO_DEPENDENCIES}}
+// Microsoft Visual C++ generated include file.
+// Used by MbDCode.rc
+
+// Next default values for new objects
+// 
+#ifdef APSTUDIO_INVOKED
+#ifndef APSTUDIO_READONLY_SYMBOLS
+#define _APS_NEXT_RESOURCE_VALUE        101
+#define _APS_NEXT_COMMAND_VALUE         40001
+#define _APS_NEXT_CONTROL_VALUE         1001
+#define _APS_NEXT_SYMED_VALUE           101
+#endif
+#endif