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

// Created on: 1998-11-17
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.


#include <FEmTool_Assembly.hxx>
#include <FEmTool_ListIteratorOfListOfVectors.hxx>
#include <FEmTool_ListOfVectors.hxx>
#include <FEmTool_ProfileMatrix.hxx>
#include <math_Matrix.hxx>
#include <Standard_DimensionError.hxx>
#include <Standard_DomainError.hxx>
#include <StdFail_NotDone.hxx>
#include <TColStd_HArray1OfReal.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();

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