[occt.git] / src / FEmTool / FEmTool_Assembly.cxx

// File:	FEmTool_Assembly.cxx
// Created:	Tue Nov 17 12:07:58 1998
// Author:	Igor FEOKTISTOV
//		<ifv@paradox.nnov.matra-dtv.fr>


#include <FEmTool_Assembly.ixx>
#include <FEmTool_ListIteratorOfListOfVectors.hxx>
#include <FEmTool_ListOfVectors.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <math_Matrix.hxx>

//----------------------------------------------------------------------------
// Purpose - to find min index of global variables and define
//----------------------------------------------------------------------------
static Standard_Integer MinIndex(const Handle(FEmTool_HAssemblyTable)& Table)
{
  Standard_Integer dim, el, nvar, Imin;
  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);

  nvarl = T->Lower();

  Imin = T->Value(nvarl);

  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {
      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	Imin = Min(Imin, T->Value(nvar));
      }
    }
  return Imin;
}

//----------------------------------------------------------------------------
// Purpose - to find max index of global variables 
//----------------------------------------------------------------------------
static Standard_Integer MaxIndex(const Handle(FEmTool_HAssemblyTable)& Table)
{
  Standard_Integer dim, el, nvar, Imax;
  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);

  nvarl = T->Lower();

  Imax = T->Value(nvarl);

  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {
      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	Imax = Max(Imax, T->Value(nvar));
      }
    }
  return Imax;
}


//=======================================================================
//function : FEmTool_Assembly
//purpose  : 
//=======================================================================
FEmTool_Assembly::FEmTool_Assembly(const TColStd_Array2OfInteger& Dependence,
				   const Handle(FEmTool_HAssemblyTable)& Table):
     myDepTable(1, Dependence.ColLength(), 1, Dependence.RowLength()),
     B(MinIndex(Table), MaxIndex(Table))
     
{
  IsSolved = Standard_False;
  myDepTable = Dependence;
  myRefTable = Table;

  TColStd_Array1OfInteger FirstIndexes(1, B.Length()); FirstIndexes.Init(B.Length());

  Standard_Integer dim, el, nvar, Imin, I0 = 1 - B.Lower(), i;

  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T;
  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {

      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      Imin = T->Value(nvarl) + I0;
     
      for(nvar = nvarl; nvar <= nvaru; nvar++) 
	Imin = Min(Imin, T->Value(nvar) + I0);
      
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	i = T->Value(nvar) + I0;
	FirstIndexes(i) = Min(FirstIndexes(i), Imin);
      }
    }


  H = new FEmTool_ProfileMatrix(FirstIndexes);

  NullifyMatrix();
  NullifyVector();
}

void FEmTool_Assembly::NullifyMatrix() 
{
  H->Init(0.);
  IsSolved = Standard_False;
}


//=======================================================================
//function : AddMatrix
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddMatrix(const Standard_Integer Element,
				 const Standard_Integer Dimension1,
				 const Standard_Integer Dimension2,
				 const math_Matrix& Mat)
{

  if(myDepTable(Dimension1, Dimension2) == 0) 
    Standard_DomainError::Raise("FEmTool_Assembly::AddMatrix");

  const TColStd_Array1OfInteger & T1 = myRefTable->Value(Dimension1,Element)->Array1();
  const TColStd_Array1OfInteger & T2 = myRefTable->Value(Dimension2,Element)->Array1();

  Standard_Integer nvarl = T1.Lower(), nvaru = Min(T1.Upper(), nvarl + Mat.RowNumber() - 1);


  Standard_Integer I, J, I0 = 1 - B.Lower(), i, ii, j,

                   i0 = Mat.LowerRow() - nvarl, j0 = Mat.LowerCol() - nvarl;

  for(i = nvarl; i <= nvaru; i++) {
    I = T1(i) + I0;
    ii = i0+i;
    for(j = nvarl; j <= i; j++) {
      J = T2(j) + I0;
      H->ChangeValue(I, J) += Mat(ii, j0+j);
    }
  }

  IsSolved = Standard_False;
}


//=======================================================================
//function : NullifyVector
//purpose  : 
//=======================================================================
void FEmTool_Assembly::NullifyVector() 
{
  B.Init(0.);
}

//=======================================================================
//function : AddVector
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddVector(const Standard_Integer Element,
				 const Standard_Integer Dimension,
				 const math_Vector& Vec) 
{
  const TColStd_Array1OfInteger & T = myRefTable->Value(Dimension,Element)->Array1();
  Standard_Integer nvarl = T.Lower(), nvaru = Min(T.Upper(), nvarl + Vec.Length() - 1),
                   i0 = Vec.Lower() - nvarl;

//  Standard_Integer I, i;
  Standard_Integer i;
  
  for(i = nvarl; i <= nvaru; i++) 
    B(T(i)) += Vec(i0 + i);
  
}


//=======================================================================
//function : Solve
//purpose  : 
//=======================================================================
Standard_Boolean FEmTool_Assembly::Solve() 
{
  IsSolved = H->Decompose();

#if DEB
  if (!IsSolved) {
    cout << "Solve Echec  H = " << endl;
    H->OutM();
  }
#endif

/* ????
  while(!IsSolved && count < 5) {
    Standard_Integer i;
    for(i = 1; i <= H->RowNumber(); i++) {
      H->ChangeValue(i,i) *= 2.;
    }
    IsSolved = H->Decompose();
    count++;
  }
*/    

// Standard_Integer count = 0;
  if(!G.IsEmpty() && IsSolved) {
    // calculating H-1 Gt
    math_Vector gi(B.Lower(), B.Upper()), qi(B.Lower(), B.Upper());

    if(GHGt.IsNull() || GHGt->RowNumber() != G.Length()) {
      TColStd_Array1OfInteger FirstIndexes(1, G.Length());
//-----------------------------------------------------------------
      TColStd_Array2OfInteger H1(1, NbGlobVar(), 1, NbGlobVar()); H1.Init(1);
      Standard_Integer i, j, k, l, BlockBeg = 1, BlockEnd;
      Standard_Boolean Block, Zero;
      for(i = 2; i <= NbGlobVar(); i++) {
	BlockEnd = i - 1;
	if(!H->IsInProfile(i, BlockEnd)) {
	  // Maybe, begin of block
	  Block = Standard_True;
	  for(j = i + 1; j <= NbGlobVar(); j++) {
	    if(H->IsInProfile(j, BlockEnd)) {
	      Block = Standard_False;
	      break;
	    }
	  }
	  if(Block) {
	    for(j = i; j <= NbGlobVar(); j++) {
	      for(k = BlockBeg; k <= BlockEnd; k++) {
		H1(j, k) = 0; H1(k, j) = 0;
	      } 
	    }
	    BlockBeg = BlockEnd + 1;
	  }
	  else i = j;
	}
      }
      
      FEmTool_ListIteratorOfListOfVectors Iter1;
      FEmTool_ListIteratorOfListOfVectors Iter2;
      for(i = 1; i <= G.Length(); i++) {
	const FEmTool_ListOfVectors& Gi = G.Value(i);
	for(j = 1; j <= i; j++) {
	  const FEmTool_ListOfVectors& Gj = G.Value(j);
	  Zero = Standard_True;
	  for(Iter1.Initialize(Gi); Iter1.More(); Iter1.Next()) {
	    const Handle(TColStd_HArray1OfReal)& a = Iter1.Value();
	    for(k = a->Lower(); k <= a->Upper(); k++) {
	      for(Iter2.Initialize(Gj); Iter2.More(); Iter2.Next()) {
		const Handle(TColStd_HArray1OfReal)& b = Iter2.Value();
		for(l = b->Lower(); l <= b->Upper(); l++) {
		  if(H1(k, l) != 0) {
		    Zero = Standard_False;
		    break;
		  }
		}
		if(!Zero) break;
	      }
	      if(!Zero) break;
	    }
	    if(!Zero) break;
	  }
	  if(!Zero) {
	    FirstIndexes(i) = j;
	    break;
	  }
	}
      }
//-----------------------------------------------------------------------
//      for(i = FirstIndexes.Lower(); i <= FirstIndexes.Upper(); i++)
//	cout << "FirstIndexes(" << i << ") = " << FirstIndexes(i) << endl;
      //      FirstIndexes.Init(1); // temporary GHGt is full matrix
      GHGt = new FEmTool_ProfileMatrix(FirstIndexes);
    }

    GHGt->Init(0.);
    Standard_Integer i, j, k;
    FEmTool_ListIteratorOfListOfVectors Iter;
    
    for(i = 1; i <= G.Length(); i++) {

      const FEmTool_ListOfVectors& L = G.Value(i);
      gi.Init(0.);
      // preparing i-th line of G (or column of Gt) 
      for(Iter.Initialize(L); Iter.More(); Iter.Next()) {

	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();

	for(j = a->Lower(); j <= a->Upper(); j++) gi(j) = a->Value(j); // gi - full line of G 

      }
                        //                                     -1 t
      H->Solve(gi, qi); // solving H*qi = gi, qi is column of H  G


      //                               -1 t
      // Calculation of product M = G H G
      // for each i all elements of i-th column of M are calculated for k >= i
      for(k = i; k <= G.Length(); k++) {

	if(GHGt->IsInProfile(k, i)) {
	  Standard_Real m = 0.; // m = M(k,i)
	
	  const FEmTool_ListOfVectors& L = G.Value(k);
	
	  for(Iter.Initialize(L); Iter.More(); Iter.Next()) {

	    const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	    for(j = a->Lower(); j <= a->Upper(); j++) m += qi(j) * a->Value(j); // scalar product of
	                                             // k-th line of G and i-th column of H-1 Gt 

	  }
	
	  GHGt->ChangeValue(k, i) = m;
	}
      }

    }


    IsSolved = GHGt->Decompose();
/*    count = 0;
    while(!IsSolved && count < 5) {
      for(i = 1; i <= GHGt->RowNumber(); i++) {
	GHGt->ChangeValue(i,i) *= 2.;
      }
      IsSolved = GHGt->Decompose();
      count++;
    }*/

  }

  return IsSolved;

}


//=======================================================================
//function : Solution
//purpose  : 
//=======================================================================
void FEmTool_Assembly::Solution(math_Vector& Solution) const
{
  if(!IsSolved) StdFail_NotDone::Raise("FEmTool_Assembly::Solution");

  if(G.IsEmpty()) H->Solve(B, Solution);
  else {

    math_Vector v1(B.Lower(), B.Upper());
    H->Solve(B, v1);
    
    math_Vector l(1, G.Length()), v2(1, G.Length());
    Standard_Integer i, j;
    FEmTool_ListIteratorOfListOfVectors Iter;
    
    for(i = 1; i <= G.Length(); i++) {
      
      const FEmTool_ListOfVectors& L = G.Value(i);
      Standard_Real m = 0.;
      
      for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
	
	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	for(j = a->Lower(); j <= a->Upper(); j++) 
	  m += v1(j) * a->Value(j); // scalar product
	// G v1
      }
      
      v2(i) = m - C.Value(i);
    }
    
    GHGt->Solve(v2, l); // Solving M*l = v2
    
    v1 = B;
    // Calculation v1 = B-Gt*l
    // v1(j) = B(j) - Gt(j,i)*l(i) = B(j) - G(i,j)*l(i)
    // 
      for(i = 1; i <= G.Length(); i++) {
	
	const FEmTool_ListOfVectors& L = G.Value(i);
	
	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
	  
	  const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	  for(j = a->Lower(); j <= a->Upper(); j++) v1(j) -= l(i) * a->Value(j);
	  
	}
	
      }
      
    H->Solve(v1, Solution);
  }
  
}

Standard_Integer FEmTool_Assembly::NbGlobVar() const
{

  return B.Length();

}

void FEmTool_Assembly::GetAssemblyTable(Handle(FEmTool_HAssemblyTable)& AssTable) const
{
  AssTable = myRefTable;
}


void FEmTool_Assembly::ResetConstraint() 
{
  G.Clear();
  C.Clear();
}

void FEmTool_Assembly::NullifyConstraint() 
{
  FEmTool_ListIteratorOfListOfVectors Iter;
  Standard_Integer i;

  for(i = 1; i <= G.Length(); i++) {
  
    C.SetValue(i, 0.);

    for(Iter.Initialize(G.Value(i)); Iter.More(); Iter.Next()) 
      Iter.Value()->Init(0.);

  }

}


//=======================================================================
//function : AddConstraint
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddConstraint(const Standard_Integer IndexofConstraint,
				     const Standard_Integer Element,
				     const Standard_Integer Dimension,
				     const math_Vector& LinearForm,
				     const Standard_Real Value) 
{
  while(G.Length() < IndexofConstraint) {
    // Add new lines in G
    FEmTool_ListOfVectors L;
    G.Append(L);
    C.Append(0.);
  }

  FEmTool_ListOfVectors& L = G.ChangeValue(IndexofConstraint);
  
  Handle(TColStd_HArray1OfInteger) Indexes = myRefTable->Value(Dimension,Element);
  Standard_Integer i, Imax = 0, Imin = NbGlobVar();

  for(i = Indexes->Lower(); i <= Indexes->Upper(); i++) {
    Imin = Min(Imin, Indexes->Value(i));
    Imax = Max(Imax, Indexes->Value(i));
  }

  Handle(TColStd_HArray1OfReal) Coeff;

  if(L.IsEmpty()) {
    Coeff = new TColStd_HArray1OfReal(Imin,Imax);
    Coeff->Init(0.);
    L.Append(Coeff);
  }
  else {
    FEmTool_ListIteratorOfListOfVectors Iter(L);
    Standard_Integer i;
    Standard_Real  s1 = 0, s2 = 0;
    Handle(TColStd_HArray1OfReal) Aux1, Aux2;
    for(i=1; Iter.More(); Iter.Next(), i++) {
      if(Imin >= Iter.Value()->Lower()) {
	s1 = i;
	Aux1 = Iter.Value();
	if(Imax <= Iter.Value()->Upper()) {
	  s2 = s1;
	  Coeff = Iter.Value();
	  break;
	}
      }

      if(Imax <= Iter.Value()->Upper()) {
	s2 = i;
	Aux2 = Iter.Value();
      }
    }

    if(s1 != s2) {
      if(s1 == 0) {
	if(Imax < Aux2->Lower()) {
	  // inserting before first segment
	  Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	  Coeff->Init(0.);
	  L.Prepend(Coeff);
	}
	else {
	  // merge new and first segment
	  Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
	  for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	  for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	  L.First() = Coeff;
	}
      }
      else if(s2 == 0) {
	if(Imin > Aux1->Upper()) {
	  // append new
	  Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	  Coeff->Init(0.);
	  L.Append(Coeff);
	}
	else {
	  // merge new and last segment
	  Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
	  for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	  for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
	  L.Last() = Coeff;
	}
      }
      else if(Imin <= Aux1->Upper() && Imax < Aux2->Lower()) {
	// merge s1 and new
	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	Iter.Value() = Coeff;
      }
      else if(Imin > Aux1->Upper() && Imax >= Aux2->Lower()) { 
	// merge new and first segment
	Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
	for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	Iter.Initialize(L);
	for(i = 1; i < s2; Iter.Next(), i++);
	Iter.Value() = Coeff;
     }
      else if(Imin > Aux1->Upper() && Imax < Aux2->Lower()) {
	// inserting new between s1 and s2
	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	Coeff->Init(0.);
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	L.InsertAfter(Coeff,Iter);
      }
      else {
	// merge s1, new, s2 and remove s2 
	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Aux2->Upper());
	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	for(i = Aux1->Upper() + 1; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	Iter.Value() = Coeff;
	Iter.Next();
	L.Remove(Iter);
      }
    }
  }

  // adding 
  Standard_Integer j = LinearForm.Lower();
  for(i = Indexes->Lower(); i <= Indexes->Upper(); i++, j++) {
    Coeff->ChangeValue(Indexes->Value(i)) += LinearForm(j);
  }
  
  C.ChangeValue(IndexofConstraint) += Value;


}
Commit	Line	Data
7fd59977	1	// File: FEmTool_Assembly.cxx
	2	// Created: Tue Nov 17 12:07:58 1998
	3	// Author: Igor FEOKTISTOV
	4	// <ifv@paradox.nnov.matra-dtv.fr>
	5
	6
	7	#include <FEmTool_Assembly.ixx>
	8	#include <FEmTool_ListIteratorOfListOfVectors.hxx>
	9	#include <FEmTool_ListOfVectors.hxx>
	10	#include <TColStd_HArray1OfReal.hxx>
	11	#include <math_Matrix.hxx>
	12
	13	//----------------------------------------------------------------------------
	14	// Purpose - to find min index of global variables and define
	15	//----------------------------------------------------------------------------
	16	static Standard_Integer MinIndex(const Handle(FEmTool_HAssemblyTable)& Table)
	17	{
	18	Standard_Integer dim, el, nvar, Imin;
	19	Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
	20	ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
	21
	22	Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);
	23
	24	nvarl = T->Lower();
	25
	26	Imin = T->Value(nvarl);
	27
	28	for(dim = diml; dim <= dimu; dim++)
	29	for(el = ell; el <= elu; el++) {
	30	T = Table->Value(dim, el);
	31	nvarl = T->Lower(); nvaru = T->Upper();
	32	for(nvar = nvarl; nvar <= nvaru; nvar++) {
	33	Imin = Min(Imin, T->Value(nvar));
	34	}
	35	}
	36	return Imin;
	37	}
	38
	39	//----------------------------------------------------------------------------
	40	// Purpose - to find max index of global variables
	41	//----------------------------------------------------------------------------
	42	static Standard_Integer MaxIndex(const Handle(FEmTool_HAssemblyTable)& Table)
	43	{
	44	Standard_Integer dim, el, nvar, Imax;
	45	Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
	46	ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
	47
	48	Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);
	49
	50	nvarl = T->Lower();
	51
	52	Imax = T->Value(nvarl);
	53
	54	for(dim = diml; dim <= dimu; dim++)
	55	for(el = ell; el <= elu; el++) {
	56	T = Table->Value(dim, el);
	57	nvarl = T->Lower(); nvaru = T->Upper();
	58	for(nvar = nvarl; nvar <= nvaru; nvar++) {
	59	Imax = Max(Imax, T->Value(nvar));
	60	}
	61	}
	62	return Imax;
	63	}
	64
65
66
67	//=======================================================================
68	//function : FEmTool_Assembly
69	//purpose :
70	//=======================================================================
71	FEmTool_Assembly::FEmTool_Assembly(const TColStd_Array2OfInteger& Dependence,
72	const Handle(FEmTool_HAssemblyTable)& Table):
73	myDepTable(1, Dependence.ColLength(), 1, Dependence.RowLength()),
74	B(MinIndex(Table), MaxIndex(Table))
75
76	{
77	IsSolved = Standard_False;
78	myDepTable = Dependence;
79	myRefTable = Table;
80
81	TColStd_Array1OfInteger FirstIndexes(1, B.Length()); FirstIndexes.Init(B.Length());
82
7fd59977	83	Standard_Integer dim, el, nvar, Imin, I0 = 1 - B.Lower(), i;
6e6cd5d9	84
7fd59977	85	Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
	86	ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
	87
	88	Handle(TColStd_HArray1OfInteger) T;
	89	for(dim = diml; dim <= dimu; dim++)
	90	for(el = ell; el <= elu; el++) {
	91
	92	T = Table->Value(dim, el);
	93	nvarl = T->Lower(); nvaru = T->Upper();
	94	Imin = T->Value(nvarl) + I0;
	95
	96	for(nvar = nvarl; nvar <= nvaru; nvar++)
	97	Imin = Min(Imin, T->Value(nvar) + I0);
	98
	99	for(nvar = nvarl; nvar <= nvaru; nvar++) {
	100	i = T->Value(nvar) + I0;
	101	FirstIndexes(i) = Min(FirstIndexes(i), Imin);
	102	}
	103	}
	104
	105
	106	H = new FEmTool_ProfileMatrix(FirstIndexes);
	107
	108	NullifyMatrix();
	109	NullifyVector();
	110	}
	111
	112	void FEmTool_Assembly::NullifyMatrix()
	113	{
	114	H->Init(0.);
	115	IsSolved = Standard_False;
	116	}
	117
	118
	119	//=======================================================================
	120	//function : AddMatrix
	121	//purpose :
	122	//=======================================================================
	123	void FEmTool_Assembly::AddMatrix(const Standard_Integer Element,
	124	const Standard_Integer Dimension1,
	125	const Standard_Integer Dimension2,
	126	const math_Matrix& Mat)
	127	{
	128
	129	if(myDepTable(Dimension1, Dimension2) == 0)
	130	Standard_DomainError::Raise("FEmTool_Assembly::AddMatrix");
	131
	132	const TColStd_Array1OfInteger & T1 = myRefTable->Value(Dimension1,Element)->Array1();
	133	const TColStd_Array1OfInteger & T2 = myRefTable->Value(Dimension2,Element)->Array1();
	134
	135	Standard_Integer nvarl = T1.Lower(), nvaru = Min(T1.Upper(), nvarl + Mat.RowNumber() - 1);
	136
6e6cd5d9	137
7fd59977	138	Standard_Integer I, J, I0 = 1 - B.Lower(), i, ii, j,
6e6cd5d9	139
7fd59977	140	i0 = Mat.LowerRow() - nvarl, j0 = Mat.LowerCol() - nvarl;
	141
	142	for(i = nvarl; i <= nvaru; i++) {
	143	I = T1(i) + I0;
	144	ii = i0+i;
	145	for(j = nvarl; j <= i; j++) {
	146	J = T2(j) + I0;
	147	H->ChangeValue(I, J) += Mat(ii, j0+j);
	148	}
	149	}
	150
	151	IsSolved = Standard_False;
	152	}
	153
	154
	155	//=======================================================================
	156	//function : NullifyVector
	157	//purpose :
	158	//=======================================================================
	159	void FEmTool_Assembly::NullifyVector()
	160	{
	161	B.Init(0.);
	162	}
	163
	164	//=======================================================================
	165	//function : AddVector
	166	//purpose :
	167	//=======================================================================
	168	void FEmTool_Assembly::AddVector(const Standard_Integer Element,
	169	const Standard_Integer Dimension,
	170	const math_Vector& Vec)
	171	{
	172	const TColStd_Array1OfInteger & T = myRefTable->Value(Dimension,Element)->Array1();
	173	Standard_Integer nvarl = T.Lower(), nvaru = Min(T.Upper(), nvarl + Vec.Length() - 1),
	174	i0 = Vec.Lower() - nvarl;
	175
	176	// Standard_Integer I, i;
	177	Standard_Integer i;
	178
	179	for(i = nvarl; i <= nvaru; i++)
	180	B(T(i)) += Vec(i0 + i);
	181
	182	}
	183
	184
	185	//=======================================================================
	186	//function : Solve
	187	//purpose :
	188	//=======================================================================
	189	Standard_Boolean FEmTool_Assembly::Solve()
	190	{
	191	IsSolved = H->Decompose();
	192
	193	#if DEB
	194	if (!IsSolved) {
	195	cout << "Solve Echec H = " << endl;
	196	H->OutM();
	197	}
	198	#endif
	199
	200	/* ????
	201	while(!IsSolved && count < 5) {
	202	Standard_Integer i;
	203	for(i = 1; i <= H->RowNumber(); i++) {
204	H->ChangeValue(i,i) *= 2.;
205	}
206	IsSolved = H->Decompose();
207	count++;
208	}
209	*/
210
211	// Standard_Integer count = 0;
212	if(!G.IsEmpty() && IsSolved) {
213	// calculating H-1 Gt
214	math_Vector gi(B.Lower(), B.Upper()), qi(B.Lower(), B.Upper());
215
216	if(GHGt.IsNull() \|\| GHGt->RowNumber() != G.Length()) {
217	TColStd_Array1OfInteger FirstIndexes(1, G.Length());
218	//-----------------------------------------------------------------
219	TColStd_Array2OfInteger H1(1, NbGlobVar(), 1, NbGlobVar()); H1.Init(1);
220	Standard_Integer i, j, k, l, BlockBeg = 1, BlockEnd;
221	Standard_Boolean Block, Zero;
222	for(i = 2; i <= NbGlobVar(); i++) {
223	BlockEnd = i - 1;
224	if(!H->IsInProfile(i, BlockEnd)) {
225	// Maybe, begin of block
226	Block = Standard_True;
227	for(j = i + 1; j <= NbGlobVar(); j++) {
228	if(H->IsInProfile(j, BlockEnd)) {
229	Block = Standard_False;
230	break;
231	}
232	}
233	if(Block) {
234	for(j = i; j <= NbGlobVar(); j++) {
235	for(k = BlockBeg; k <= BlockEnd; k++) {
236	H1(j, k) = 0; H1(k, j) = 0;
237	}
238	}
239	BlockBeg = BlockEnd + 1;
240	}
241	else i = j;
242	}
243	}
244
245	FEmTool_ListIteratorOfListOfVectors Iter1;
246	FEmTool_ListIteratorOfListOfVectors Iter2;
247	for(i = 1; i <= G.Length(); i++) {
248	const FEmTool_ListOfVectors& Gi = G.Value(i);
249	for(j = 1; j <= i; j++) {
250	const FEmTool_ListOfVectors& Gj = G.Value(j);
251	Zero = Standard_True;
252	for(Iter1.Initialize(Gi); Iter1.More(); Iter1.Next()) {
253	const Handle(TColStd_HArray1OfReal)& a = Iter1.Value();
254	for(k = a->Lower(); k <= a->Upper(); k++) {
255	for(Iter2.Initialize(Gj); Iter2.More(); Iter2.Next()) {
256	const Handle(TColStd_HArray1OfReal)& b = Iter2.Value();
257	for(l = b->Lower(); l <= b->Upper(); l++) {
258	if(H1(k, l) != 0) {
259	Zero = Standard_False;
260	break;
261	}
262	}
263	if(!Zero) break;
264	}
265	if(!Zero) break;
266	}
267	if(!Zero) break;
268	}
269	if(!Zero) {
270	FirstIndexes(i) = j;
271	break;
272	}
273	}
274	}
275	//-----------------------------------------------------------------------
276	// for(i = FirstIndexes.Lower(); i <= FirstIndexes.Upper(); i++)
277	// cout << "FirstIndexes(" << i << ") = " << FirstIndexes(i) << endl;
278	// FirstIndexes.Init(1); // temporary GHGt is full matrix
279	GHGt = new FEmTool_ProfileMatrix(FirstIndexes);
280	}
281
282	GHGt->Init(0.);
283	Standard_Integer i, j, k;
284	FEmTool_ListIteratorOfListOfVectors Iter;
285
286	for(i = 1; i <= G.Length(); i++) {
287
288	const FEmTool_ListOfVectors& L = G.Value(i);
289	gi.Init(0.);
290	// preparing i-th line of G (or column of Gt)
291	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
292
293	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
294
295	for(j = a->Lower(); j <= a->Upper(); j++) gi(j) = a->Value(j); // gi - full line of G
296
297	}
298	// -1 t
299	H->Solve(gi, qi); // solving H*qi = gi, qi is column of H G
300
301
302	// -1 t
303	// Calculation of product M = G H G
304	// for each i all elements of i-th column of M are calculated for k >= i
305	for(k = i; k <= G.Length(); k++) {
306
307	if(GHGt->IsInProfile(k, i)) {
308	Standard_Real m = 0.; // m = M(k,i)
309
310	const FEmTool_ListOfVectors& L = G.Value(k);
311
312	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
313
314	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
315	for(j = a->Lower(); j <= a->Upper(); j++) m += qi(j) * a->Value(j); // scalar product of
316	// k-th line of G and i-th column of H-1 Gt
317
318	}
319
320	GHGt->ChangeValue(k, i) = m;
321	}
322	}
323
324	}
325
326
327	IsSolved = GHGt->Decompose();
328	/* count = 0;
329	while(!IsSolved && count < 5) {
330	for(i = 1; i <= GHGt->RowNumber(); i++) {
331	GHGt->ChangeValue(i,i) *= 2.;
332	}
333	IsSolved = GHGt->Decompose();
334	count++;
335	}*/
336
337	}
338
339	return IsSolved;
340
341	}
342
343
344	//=======================================================================
345	//function : Solution
346	//purpose :
347	//=======================================================================
348	void FEmTool_Assembly::Solution(math_Vector& Solution) const
349	{
350	if(!IsSolved) StdFail_NotDone::Raise("FEmTool_Assembly::Solution");
351
352	if(G.IsEmpty()) H->Solve(B, Solution);
353	else {
354
355	math_Vector v1(B.Lower(), B.Upper());
356	H->Solve(B, v1);
357
358	math_Vector l(1, G.Length()), v2(1, G.Length());
359	Standard_Integer i, j;
360	FEmTool_ListIteratorOfListOfVectors Iter;
361
362	for(i = 1; i <= G.Length(); i++) {
363
364	const FEmTool_ListOfVectors& L = G.Value(i);
365	Standard_Real m = 0.;
366
367	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
368
369	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
370	for(j = a->Lower(); j <= a->Upper(); j++)
371	m += v1(j) * a->Value(j); // scalar product
372	// G v1
373	}
374
375	v2(i) = m - C.Value(i);
376	}
377
378	GHGt->Solve(v2, l); // Solving M*l = v2
379
380	v1 = B;
381	// Calculation v1 = B-Gt*l
382	// v1(j) = B(j) - Gt(j,i)l(i) = B(j) - G(i,j)l(i)
383	//
384	for(i = 1; i <= G.Length(); i++) {
385
386	const FEmTool_ListOfVectors& L = G.Value(i);
387
388	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
389
390	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
391	for(j = a->Lower(); j <= a->Upper(); j++) v1(j) -= l(i) * a->Value(j);
392
393	}
394
395	}
396
397	H->Solve(v1, Solution);
398	}
399
400	}
401
402	Standard_Integer FEmTool_Assembly::NbGlobVar() const
403	{
404
405	return B.Length();
406
407	}
408
409	void FEmTool_Assembly::GetAssemblyTable(Handle(FEmTool_HAssemblyTable)& AssTable) const
410	{
411	AssTable = myRefTable;
412	}
413
414
415	void FEmTool_Assembly::ResetConstraint()
416	{
417	G.Clear();
418	C.Clear();
419	}
420
421	void FEmTool_Assembly::NullifyConstraint()
422	{
423	FEmTool_ListIteratorOfListOfVectors Iter;
424	Standard_Integer i;
425
426	for(i = 1; i <= G.Length(); i++) {
427
428	C.SetValue(i, 0.);
429
430	for(Iter.Initialize(G.Value(i)); Iter.More(); Iter.Next())
431	Iter.Value()->Init(0.);
432
433	}
434
435	}
436
437
438	//=======================================================================
439	//function : AddConstraint
440	//purpose :
441	//=======================================================================
442	void FEmTool_Assembly::AddConstraint(const Standard_Integer IndexofConstraint,
443	const Standard_Integer Element,
444	const Standard_Integer Dimension,
445	const math_Vector& LinearForm,
446	const Standard_Real Value)
447	{
448	while(G.Length() < IndexofConstraint) {
449	// Add new lines in G
450	FEmTool_ListOfVectors L;
451	G.Append(L);
452	C.Append(0.);
453	}
454
455	FEmTool_ListOfVectors& L = G.ChangeValue(IndexofConstraint);
456
457	Handle(TColStd_HArray1OfInteger) Indexes = myRefTable->Value(Dimension,Element);
458	Standard_Integer i, Imax = 0, Imin = NbGlobVar();
459
460	for(i = Indexes->Lower(); i <= Indexes->Upper(); i++) {
461	Imin = Min(Imin, Indexes->Value(i));
462	Imax = Max(Imax, Indexes->Value(i));
463	}
464
465	Handle(TColStd_HArray1OfReal) Coeff;
466
467	if(L.IsEmpty()) {
468	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
469	Coeff->Init(0.);
470	L.Append(Coeff);
471	}
472	else {
473	FEmTool_ListIteratorOfListOfVectors Iter(L);
474	Standard_Integer i;
475	Standard_Real s1 = 0, s2 = 0;
476	Handle(TColStd_HArray1OfReal) Aux1, Aux2;
477	for(i=1; Iter.More(); Iter.Next(), i++) {
478	if(Imin >= Iter.Value()->Lower()) {
479	s1 = i;
480	Aux1 = Iter.Value();
481	if(Imax <= Iter.Value()->Upper()) {
482	s2 = s1;
483	Coeff = Iter.Value();
484	break;
485	}
486	}
487
488	if(Imax <= Iter.Value()->Upper()) {
489	s2 = i;
490	Aux2 = Iter.Value();
491	}
492	}
493
494	if(s1 != s2) {
495	if(s1 == 0) {
496	if(Imax < Aux2->Lower()) {
497	// inserting before first segment
498	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
499	Coeff->Init(0.);
500	L.Prepend(Coeff);
501	}
502	else {
503	// merge new and first segment
504	Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
505	for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
506	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
507	L.First() = Coeff;
508	}
509	}
510	else if(s2 == 0) {
511	if(Imin > Aux1->Upper()) {
512	// append new
513	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
514	Coeff->Init(0.);
515	L.Append(Coeff);
516	}
517	else {
518	// merge new and last segment
519	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
520	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
521	for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
522	L.Last() = Coeff;
523	}
524	}
525	else if(Imin <= Aux1->Upper() && Imax < Aux2->Lower()) {
526	// merge s1 and new
527	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
528	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
529	for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
530	Iter.Initialize(L);
531	for(i = 1; i < s1; Iter.Next(), i++);
532	Iter.Value() = Coeff;
533	}
534	else if(Imin > Aux1->Upper() && Imax >= Aux2->Lower()) {
535	// merge new and first segment
536	Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
537	for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
538	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
539	Iter.Initialize(L);
540	for(i = 1; i < s2; Iter.Next(), i++);
541	Iter.Value() = Coeff;
542	}
543	else if(Imin > Aux1->Upper() && Imax < Aux2->Lower()) {
544	// inserting new between s1 and s2
545	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
546	Coeff->Init(0.);
547	Iter.Initialize(L);
548	for(i = 1; i < s1; Iter.Next(), i++);
549	L.InsertAfter(Coeff,Iter);
550	}
551	else {
552	// merge s1, new, s2 and remove s2
553	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Aux2->Upper());
554	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
555	for(i = Aux1->Upper() + 1; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
556	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
557	Iter.Initialize(L);
558	for(i = 1; i < s1; Iter.Next(), i++);
559	Iter.Value() = Coeff;
560	Iter.Next();
561	L.Remove(Iter);
562	}
563	}
564	}
565
566	// adding
567	Standard_Integer j = LinearForm.Lower();
568	for(i = Indexes->Lower(); i <= Indexes->Upper(); i++, j++) {
569	Coeff->ChangeValue(Indexes->Value(i)) += LinearForm(j);
570	}
571
572	C.ChangeValue(IndexofConstraint) += Value;
573
574
575
576	}
577