[occt.git] / src / AppParCurves / AppParCurves_Variational_8.gxx

// Created on: 1999-02-15
// Created by: Igor FEOKTISTOV
// Copyright (c) 1999-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.


#include <PLib_Base.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <math_Vector.hxx>

void AppParCurves_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
						     const TColStd_Array1OfReal& Parameters,
						     const Standard_Real CBLONG,
						     FEmTool_Assembly& A ) const
{

  Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
                   NbElm = Curve->NbElements(),
                   NbDim = Curve->Dimension();

  TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg);
  math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg), 
              V1((Standard_Real*)&G1(0), 0, MxDeg), 
              V2((Standard_Real*)&G2(0), 0, MxDeg);

  Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1, 
                   NTang3d, NTang2d,  
                   Point, TypOfConstr, 
                   p0 = Parameters.Lower() - myFirstPoint,
                   curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim,
                   jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d;
  Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
  
//  Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints;
//  Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints;
  Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints;
  Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints;
  NBeg2d = Ng3d * myNbP3d;
//  NgPC1 = NbConstr + myNbCurvPoints;
  NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints;
  NPass = 0; 
  NTang3d = 3 * NgPC1;
  NTang2d = 2 * NgPC1;

  TColStd_Array1OfReal& Intervals = Curve->Knots();

  Standard_Real t, R1, R2;

  Handle(PLib_Base) myBase = Curve->Base();
  Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));
  Standard_Integer Order = myHermitJacobi->NivConstr() + 1;

  Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1;

  A.NullifyConstraint();

  ipnt = -1;
  ityp = 0;
  for(i = 1; i <= NbConstr; i++) {

    ipnt += 2; ityp += 2;

    Point = myTypConstraints->Value(ipnt);
    TypOfConstr = myTypConstraints->Value(ityp);

    t = Parameters(p0 + Point);

    for(el = curel; el <= NbElm; ) {
      if( t <= Intervals(++el) ) {
	curel = el - 1;
	break;
      }
    }
    

    UFirst = Intervals(curel); ULast = Intervals(curel + 1);
    coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.;

    t = (t - c0) / coeff;
    
    if(TypOfConstr == 0) {
      myBase->D0(t, G0);
      for(k = 1; k < Order; k++) {
	mfact = Pow(coeff, k);
	G0(k) *= mfact;
	G0(k + Order) *= mfact;
      }
    }
    else if(TypOfConstr == 1) {
      myBase->D1(t, G0, G1);
      for(k = 1; k < Order; k++) {
	mfact = Pow(coeff, k);
	G0(k) *= mfact;
	G0(k + Order) *= mfact;
	G1(k) *= mfact;
	G1(k + Order) *= mfact;
      }
      mfact = 1./coeff;
      for(k = 0; k <= MxDeg; k++) {
	G1(k) *= mfact;
      }
    }
    else {
      myBase->D2(t, G0, G1, G2);
      for(k = 1; k < Order; k++) {
	mfact = Pow(coeff, k);
	G0(k) *= mfact;
	G0(k + Order) *= mfact;
	G1(k) *= mfact;
	G1(k + Order) *= mfact;
	G2(k) *= mfact;
	G2(k + Order) *= mfact;
      }
      mfact = 1. / coeff;
      mfact1 = mfact / coeff;
      for(k = 0; k <= MxDeg; k++) {
	G1(k) *= mfact;
	G2(k) *= mfact1;
      }
    }

    NPass++;

    j = NbDim * (Point - myFirstPoint);
    Standard_Integer n0 = NPass;
    curdim = 0;
    for(pnt = 1; pnt <= myNbP3d; pnt++) {
      IndexOfConstraint = n0;
      for(k = 1; k <= 3; k++) {
	curdim++;
	A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
	IndexOfConstraint += NgPC1;
      }
      j +=3;
      n0 += Ng3d;
    }

    n0 = NPass + NBeg2d;
    for(pnt = 1; pnt <= myNbP2d; pnt++) {
      IndexOfConstraint = n0;
      for(k = 1; k <= 2; k++) {
	curdim++;
	A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
	IndexOfConstraint += NgPC1;
      }
      j +=2;
      n0 += Ng2d;
    }

/*    if(TypOfConstr == 1) {

      IndexOfConstraint = NTang3d + 1;
      jt = Ntheta * (i - 1);
      for(pnt = 1; pnt <= myNbP3d; pnt++) {
	for(k = 1; k <= 3; k++) {
	  A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) *  V1, 0.);
	  A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) *  V1, 0.);
	}
	IndexOfConstraint += Ng3d;
	jt += 6;
      }

      IndexOfConstraint = NBeg2d + NTang2d + 1;
      for(pnt = 1; pnt <= myNbP2d; pnt++) {
	for(k = 1; k <= 2; k++) {
	  A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
	}
	IndexOfConstraint += Ng2d;
	jt += 2;
      }

      NTang3d += 2;
      NTang2d += 1;
    } */
    if(TypOfConstr == 1) {

      NPass++;
      n0 = NPass;
      j = 2 * NbDim * (i - 1);
      curdim = 0;
      for(pnt = 1; pnt <= myNbP3d; pnt++) {
	IndexOfConstraint = n0;
	for(k = 1; k <= 3; k++) {
	  curdim++;
	  A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
	  IndexOfConstraint += NgPC1;
	}
	n0 += Ng3d;
	j += 6;
      }

      n0 = NPass + NBeg2d;
      for(pnt = 1; pnt <= myNbP2d; pnt++) {
	IndexOfConstraint = n0;
	for(k = 1; k <= 2; k++) {
	  curdim++;
	  A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
	  IndexOfConstraint += NgPC1;
	}
	n0 += Ng2d;
	j += 4;
      }
    }
    if(TypOfConstr == 2) {

      NPass++;
      n0 = NPass;
      j = 2 * NbDim * (i - 1);
      curdim = 0;
      for(pnt = 1; pnt <= myNbP3d; pnt++) {
	IndexOfConstraint = n0;
	for(k = 1; k <= 3; k++) {
	  curdim++;
	  A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
	  IndexOfConstraint += NgPC1;
	}
	n0 += Ng3d;
	j += 6;
      }

      n0 = NPass + NBeg2d;
      for(pnt = 1; pnt <= myNbP2d; pnt++) {
	IndexOfConstraint = n0;
	for(k = 1; k <= 2; k++) {
	  curdim++;
	  A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
	  IndexOfConstraint += NgPC1;
	}
	n0 += Ng2d;
	j += 4;
      }

      j = 2 * NbDim * (i - 1) + 3;
      jt = Ntheta * (i - 1);
      IndexOfConstraint = NTang3d + 1;
      curdim = 0;
      for(pnt = 1; pnt <= myNbP3d; pnt++) {
	R1 = 0.; R2 = 0.;
	for(k = 1; k <= 3; k++) { 
	  R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
	  R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k);
	}
	R1 *= CBLONG * CBLONG;
	R2 *= CBLONG * CBLONG;
	for(k = 1; k <= 3; k++) {
	  curdim++;
	  if(k > 1) R1 = R2 = 0.;
	  A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
	  A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2);
	}
	IndexOfConstraint += Ng3d;
	j += 6;
	jt += 6;
      }

      j--;
      IndexOfConstraint = NBeg2d + NTang2d + 1;
      for(pnt = 1; pnt <= myNbP2d; pnt++) {
	R1 = 0.;
	for(k = 1; k <= 2; k++) { 
	  R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
	}
	R1 *= CBLONG * CBLONG;
	for(k = 1; k <= 2; k++) {
	  curdim++;
	  if(k > 1) R1 = 0.;
	  A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
	}
	IndexOfConstraint += Ng2d;
	j += 4;
	jt += 2;
      }

      NTang3d += 2;
      NTang2d += 1;
    }

  }

}
Commit	Line	Data
b311480e	1	// Created on: 1999-02-15
	2	// Created by: Igor FEOKTISTOV
	3	// Copyright (c) 1999-1999 Matra Datavision
	4	// Copyright (c) 1999-2012 OPEN CASCADE SAS
	5	//
	6	// The content of this file is subject to the Open CASCADE Technology Public
	7	// License Version 6.5 (the "License"). You may not use the content of this file
	8	// except in compliance with the License. Please obtain a copy of the License
	9	// at http://www.opencascade.org and read it completely before using this file.
	10	//
	11	// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
	12	// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
	13	//
	14	// The Original Code and all software distributed under the License is
	15	// distributed on an "AS IS" basis, without warranty of any kind, and the
	16	// Initial Developer hereby disclaims all such warranties, including without
	17	// limitation, any warranties of merchantability, fitness for a particular
	18	// purpose or non-infringement. Please see the License for the specific terms
	19	// and conditions governing the rights and limitations under the License.
	20
7fd59977	21
	22	#include <PLib_Base.hxx>
	23	#include <PLib_JacobiPolynomial.hxx>
	24	#include <PLib_HermitJacobi.hxx>
	25	#include <TColStd_Array1OfReal.hxx>
	26	#include <math_Vector.hxx>
	27
	28	void AppParCurves_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
	29	const TColStd_Array1OfReal& Parameters,
	30	const Standard_Real CBLONG,
	31	FEmTool_Assembly& A ) const
	32	{
	33
	34	Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
	35	NbElm = Curve->NbElements(),
	36	NbDim = Curve->Dimension();
	37
	38	TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg);
	39	math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg),
	40	V1((Standard_Real*)&G1(0), 0, MxDeg),
	41	V2((Standard_Real*)&G2(0), 0, MxDeg);
	42
	43	Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1,
	44	NTang3d, NTang2d,
	45	Point, TypOfConstr,
	46	p0 = Parameters.Lower() - myFirstPoint,
	47	curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim,
	48	jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d;
	49	Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
	50
	51	// Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints;
	52	// Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints;
	53	Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints;
	54	Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints;
	55	NBeg2d = Ng3d * myNbP3d;
	56	// NgPC1 = NbConstr + myNbCurvPoints;
	57	NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints;
	58	NPass = 0;
	59	NTang3d = 3 * NgPC1;
	60	NTang2d = 2 * NgPC1;
	61
	62	TColStd_Array1OfReal& Intervals = Curve->Knots();
	63
	64	Standard_Real t, R1, R2;
	65
	66	Handle(PLib_Base) myBase = Curve->Base();
	67	Handle(PLib_HermitJacobi) myHermitJacobi = (((Handle(PLib_HermitJacobi))&myBase));
	68	Standard_Integer Order = myHermitJacobi->NivConstr() + 1;
	69
	70	Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1;
	71
	72	A.NullifyConstraint();
	73
	74	ipnt = -1;
	75	ityp = 0;
	76	for(i = 1; i <= NbConstr; i++) {
	77
	78	ipnt += 2; ityp += 2;
	79
	80	Point = myTypConstraints->Value(ipnt);
	81	TypOfConstr = myTypConstraints->Value(ityp);
	82
	83	t = Parameters(p0 + Point);
	84
85	for(el = curel; el <= NbElm; ) {
86	if( t <= Intervals(++el) ) {
87	curel = el - 1;
88	break;
89	}
90	}
91
92
93	UFirst = Intervals(curel); ULast = Intervals(curel + 1);
94	coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.;
95
96	t = (t - c0) / coeff;
97
98	if(TypOfConstr == 0) {
99	myBase->D0(t, G0);
100	for(k = 1; k < Order; k++) {
101	mfact = Pow(coeff, k);
102	G0(k) *= mfact;
103	G0(k + Order) *= mfact;
104	}
105	}
106	else if(TypOfConstr == 1) {
107	myBase->D1(t, G0, G1);
108	for(k = 1; k < Order; k++) {
109	mfact = Pow(coeff, k);
110	G0(k) *= mfact;
111	G0(k + Order) *= mfact;
112	G1(k) *= mfact;
113	G1(k + Order) *= mfact;
114	}
115	mfact = 1./coeff;
116	for(k = 0; k <= MxDeg; k++) {
117	G1(k) *= mfact;
118	}
119	}
120	else {
121	myBase->D2(t, G0, G1, G2);
122	for(k = 1; k < Order; k++) {
123	mfact = Pow(coeff, k);
124	G0(k) *= mfact;
125	G0(k + Order) *= mfact;
126	G1(k) *= mfact;
127	G1(k + Order) *= mfact;
128	G2(k) *= mfact;
129	G2(k + Order) *= mfact;
130	}
131	mfact = 1. / coeff;
132	mfact1 = mfact / coeff;
133	for(k = 0; k <= MxDeg; k++) {
134	G1(k) *= mfact;
135	G2(k) *= mfact1;
136	}
137	}
138
139	NPass++;
140
141	j = NbDim * (Point - myFirstPoint);
142	Standard_Integer n0 = NPass;
143	curdim = 0;
144	for(pnt = 1; pnt <= myNbP3d; pnt++) {
145	IndexOfConstraint = n0;
146	for(k = 1; k <= 3; k++) {
147	curdim++;
148	A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
149	IndexOfConstraint += NgPC1;
150	}
151	j +=3;
152	n0 += Ng3d;
153	}
154
155	n0 = NPass + NBeg2d;
156	for(pnt = 1; pnt <= myNbP2d; pnt++) {
157	IndexOfConstraint = n0;
158	for(k = 1; k <= 2; k++) {
159	curdim++;
160	A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
161	IndexOfConstraint += NgPC1;
162	}
163	j +=2;
164	n0 += Ng2d;
165	}
166
167	/* if(TypOfConstr == 1) {
168
169	IndexOfConstraint = NTang3d + 1;
170	jt = Ntheta * (i - 1);
171	for(pnt = 1; pnt <= myNbP3d; pnt++) {
172	for(k = 1; k <= 3; k++) {
173	A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
174	A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) * V1, 0.);
175	}
176	IndexOfConstraint += Ng3d;
177	jt += 6;
178	}
179
180	IndexOfConstraint = NBeg2d + NTang2d + 1;
181	for(pnt = 1; pnt <= myNbP2d; pnt++) {
182	for(k = 1; k <= 2; k++) {
183	A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
184	}
185	IndexOfConstraint += Ng2d;
186	jt += 2;
187	}
188
189	NTang3d += 2;
190	NTang2d += 1;
191	} */
192	if(TypOfConstr == 1) {
193
194	NPass++;
195	n0 = NPass;
196	j = 2 * NbDim * (i - 1);
197	curdim = 0;
198	for(pnt = 1; pnt <= myNbP3d; pnt++) {
199	IndexOfConstraint = n0;
200	for(k = 1; k <= 3; k++) {
201	curdim++;
202	A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
203	IndexOfConstraint += NgPC1;
204	}
205	n0 += Ng3d;
206	j += 6;
207	}
208
209	n0 = NPass + NBeg2d;
210	for(pnt = 1; pnt <= myNbP2d; pnt++) {
211	IndexOfConstraint = n0;
212	for(k = 1; k <= 2; k++) {
213	curdim++;
214	A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
215	IndexOfConstraint += NgPC1;
216	}
217	n0 += Ng2d;
218	j += 4;
219	}
220	}
221	if(TypOfConstr == 2) {
222
223	NPass++;
224	n0 = NPass;
225	j = 2 * NbDim * (i - 1);
226	curdim = 0;
227	for(pnt = 1; pnt <= myNbP3d; pnt++) {
228	IndexOfConstraint = n0;
229	for(k = 1; k <= 3; k++) {
230	curdim++;
231	A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
232	IndexOfConstraint += NgPC1;
233	}
234	n0 += Ng3d;
235	j += 6;
236	}
237
238	n0 = NPass + NBeg2d;
239	for(pnt = 1; pnt <= myNbP2d; pnt++) {
240	IndexOfConstraint = n0;
241	for(k = 1; k <= 2; k++) {
242	curdim++;
243	A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
244	IndexOfConstraint += NgPC1;
245	}
246	n0 += Ng2d;
247	j += 4;
248	}
249
250	j = 2 * NbDim * (i - 1) + 3;
251	jt = Ntheta * (i - 1);
252	IndexOfConstraint = NTang3d + 1;
253	curdim = 0;
254	for(pnt = 1; pnt <= myNbP3d; pnt++) {
255	R1 = 0.; R2 = 0.;
256	for(k = 1; k <= 3; k++) {
257	R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
258	R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k);
259	}
260	R1 = CBLONG CBLONG;
261	R2 = CBLONG CBLONG;
262	for(k = 1; k <= 3; k++) {
263	curdim++;
264	if(k > 1) R1 = R2 = 0.;
265	A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
266	A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2);
267	}
268	IndexOfConstraint += Ng3d;
269	j += 6;
270	jt += 6;
271	}
272
273	j--;
274	IndexOfConstraint = NBeg2d + NTang2d + 1;
275	for(pnt = 1; pnt <= myNbP2d; pnt++) {
276	R1 = 0.;
277	for(k = 1; k <= 2; k++) {
278	R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
279	}
280	R1 = CBLONG CBLONG;
281	for(k = 1; k <= 2; k++) {
282	curdim++;
283	if(k > 1) R1 = 0.;
284	A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
285	}
286	IndexOfConstraint += Ng2d;
287	j += 4;
288	jt += 2;
289	}
290
291	NTang3d += 2;
292	NTang2d += 1;
293	}
294
295	}
296
297	}