1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/struct/pzbeam2.cc,v 1.30 2017/01/12 14:46:43 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31
32 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
33
34 #include "mynewmem.h"
35 #include "pzbeam2.h"
36
37 /* PiezoActuatorBeam2 - begin */
38
39 /* Funzioni di calcolo delle matrici */
40 void
AssStiffnessMat(FullSubMatrixHandler & WMA,FullSubMatrixHandler & WMB,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)41 PiezoActuatorBeam2::AssStiffnessMat(FullSubMatrixHandler& WMA,
42 FullSubMatrixHandler& WMB,
43 doublereal dCoef,
44 const VectorHandler& XCurr,
45 const VectorHandler& XPrimeCurr)
46 {
47 DEBUGCOUT("Entering PiezoActuatorBeam2::AssStiffnessMat()" << std::endl);
48
49 Beam2::AssStiffnessMat(WMA, WMB, dCoef, XCurr, XPrimeCurr);
50
51 Mat3xN tmp1(iNumElec);
52 Mat3xN tmp2(iNumElec);
53
54 tmp1.LeftMult(R*dCoef, PiezoMat[STRAIN]);
55 WMA.Sub(1, 13, tmp1);
56 WMA.Add(7, 13, tmp1);
57 tmp2.LeftMult(Mat3x3(MatCross, p - pNode[NODE1]->GetXCurr()), tmp1);
58 WMA.Add(4, 13, tmp2);
59 tmp2.LeftMult(Mat3x3(MatCross, p - pNode[NODE2]->GetXCurr()), tmp1);
60 WMA.Sub(10, 13, tmp2);
61
62 tmp1.LeftMult(R*dCoef, PiezoMat[CURVAT]);
63 WMA.Sub(4, 13, tmp1);
64 WMA.Add(10, 13, tmp1);
65 }
66
67
68 void
AssStiffnessVec(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)69 PiezoActuatorBeam2::AssStiffnessVec(SubVectorHandler& WorkVec,
70 doublereal dCoef,
71 const VectorHandler& XCurr,
72 const VectorHandler& XPrimeCurr)
73 {
74 DEBUGCOUT("Entering PiezoActuatorBeam2::AssStiffnessVec()" << std::endl);
75
76 /*
77 * Riceve il vettore gia' dimensionato e con gli indici a posto
78 * per scrivere il residuo delle equazioni di equilibrio dei tre nodi
79 */
80
81 /*
82 * Per la trattazione teorica, il riferimento e' il file ul-travi.tex
83 * (ora e' superato)
84 */
85
86 if (bFirstRes) {
87 /*
88 * AfterPredict ha gia' calcolato tutto;
89 * bFirstRes viene resettato poi da Beam2::AssStiffnessVec
90 */
91
92 } else {
93 for (integer iCnt = 1; iCnt <= iNumElec; iCnt++) {
94 V.Put(iCnt, pvElecDofs[iCnt-1]->dGetX());
95 }
96 }
97
98 Beam2::AssStiffnessVec(WorkVec, dCoef, XCurr, XPrimeCurr);
99 }
100
101
102 void
AddInternalForces(Vec6 & AzLoc)103 PiezoActuatorBeam2::AddInternalForces(Vec6& AzLoc)
104 {
105 AzLoc += Vec6(PiezoMat[STRAIN]*V, PiezoMat[CURVAT]*V);
106 }
107
108
109 /* Costruttore normale */
PiezoActuatorBeam2(unsigned int uL,const StructNode * pN1,const StructNode * pN2,const Vec3 & F1,const Vec3 & F2,const Mat3x3 & R1,const Mat3x3 & R2,const Mat3x3 & r,const ConstitutiveLaw6D * pd,int iEl,const ScalarDifferentialNode ** pEDof,const Mat3xN & Te,const Mat3xN & Tk,OrientationDescription ood,flag fOut)110 PiezoActuatorBeam2::PiezoActuatorBeam2(unsigned int uL,
111 const StructNode* pN1, const StructNode* pN2,
112 const Vec3& F1, const Vec3& F2,
113 const Mat3x3& R1, const Mat3x3& R2,
114 const Mat3x3& r,
115 const ConstitutiveLaw6D* pd,
116 int iEl,
117 const ScalarDifferentialNode **pEDof,
118 const Mat3xN& Te, const Mat3xN& Tk,
119 OrientationDescription ood,
120 flag fOut)
121 : Elem(uL, fOut),
122 Beam2(uL, pN1, pN2, F1, F2, R1, R2, r, pd, ood, fOut),
123 iNumElec(iEl), pvElecDofs(pEDof), V(iEl)
124 {
125 ASSERT(iNumElec > 0);
126 ASSERT(pvElecDofs != NULL);
127
128 #ifdef DEBUG
129 for (int i = iNumElec; i-- > 0; ) {
130 ASSERT(pvElecDofs[i] != NULL);
131 }
132 #endif /* DEBUG */
133
134 PiezoMat[STRAIN].Copy(Te);
135 PiezoMat[CURVAT].Copy(Tk);
136 }
137
138
139 /* Distruttore banale */
~PiezoActuatorBeam2(void)140 PiezoActuatorBeam2::~PiezoActuatorBeam2(void)
141 {
142 ASSERT(pvElecDofs != NULL);
143 SAFEDELETEARR(pvElecDofs);
144 }
145
146
147 /* Contributo al file di restart */
148 std::ostream&
Restart(std::ostream & out) const149 PiezoActuatorBeam2::Restart(std::ostream& out) const
150 {
151 Restart_(out);
152 return out << "/* piezoelectric actuator NOT IMPLEMENTED YET */ ;"
153 << std::endl;
154 }
155
156
157 /*
158 * Dimensioni del workspace; sono 36 righe perche' se genera anche le
159 * forze d'inerzia consistenti deve avere accesso alle righe di definizione
160 * della quantita' di moto
161 */
162 void
WorkSpaceDim(integer * piNumRows,integer * piNumCols) const163 PiezoActuatorBeam2::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
164 {
165 *piNumRows = 12;
166 *piNumCols = 18+iNumElec;
167 }
168
169
170 /* Settings iniziali, prima della prima soluzione */
171 void
SetValue(DataManager * pDM,VectorHandler & X,VectorHandler & XP,SimulationEntity::Hints * ph)172 PiezoActuatorBeam2::SetValue(DataManager *pDM,
173 VectorHandler& X, VectorHandler& XP,
174 SimulationEntity::Hints *ph)
175 {
176 /* se proprio non serve, si puo' eliminare */
177 Beam2::SetValue(pDM, X, XP, ph);
178 }
179
180
181 /* Prepara i parametri di riferimento dopo la predizione */
182 void
AfterPredict(VectorHandler & X,VectorHandler & XP)183 PiezoActuatorBeam2::AfterPredict(VectorHandler& X, VectorHandler& XP)
184 {
185 /*
186 * Calcola le deformazioni, aggiorna il legame costitutivo
187 * e crea la FDE
188 */
189 for (integer iCnt = 1; iCnt <= iNumElec; iCnt++) {
190 V.Put(iCnt, pvElecDofs[iCnt-1]->dGetX());
191 }
192
193 Beam2::AfterPredict(X, XP);
194 }
195
196
197 /* assemblaggio jacobiano */
198 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)199 PiezoActuatorBeam2::AssJac(VariableSubMatrixHandler& WorkMat,
200 doublereal dCoef,
201 const VectorHandler& XCurr,
202 const VectorHandler& XPrimeCurr)
203 {
204 DEBUGCOUT("Entering PiezoActuatorBeam2::AssJac();"
205 " will result in call to AssStiffnessMat()" << std::endl);
206
207 integer iNode1FirstMomIndex = pNode[NODE1]->iGetFirstMomentumIndex();
208 integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
209 integer iNode2FirstMomIndex = pNode[NODE2]->iGetFirstMomentumIndex();
210 integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
211
212 FullSubMatrixHandler& WM = WorkMat.SetFull();
213
214 /* Dimensiona la matrice, la azzera e pone gli indici corretti */
215 WM.ResizeReset(12, 12+iNumElec);
216
217 for (int iCnt = 1; iCnt <= 6; iCnt++) {
218 WM.PutRowIndex(iCnt, iNode1FirstMomIndex+iCnt);
219 WM.PutColIndex(iCnt, iNode1FirstPosIndex+iCnt);
220 WM.PutRowIndex(6+iCnt, iNode2FirstMomIndex+iCnt);
221 WM.PutColIndex(6+iCnt, iNode2FirstPosIndex+iCnt);
222 }
223
224 for (int iCnt = 1; iCnt <= iNumElec; iCnt++) {
225 WM.PutColIndex(12+iCnt,
226 pvElecDofs[iCnt-1]->iGetFirstColIndex()+1);
227 }
228
229 AssStiffnessMat(WM, WM, dCoef, XCurr, XPrimeCurr);
230
231 return WorkMat;
232 }
233
234
235 /* Contributo allo jacobiano durante l'assemblaggio iniziale */
236 VariableSubMatrixHandler&
InitialAssJac(VariableSubMatrixHandler & WorkMat,const VectorHandler & XCurr)237 PiezoActuatorBeam2::InitialAssJac(VariableSubMatrixHandler& WorkMat,
238 const VectorHandler& XCurr)
239 {
240 return Beam2::InitialAssJac(WorkMat, XCurr);
241 }
242
243
244 /* Contributo al residuo durante l'assemblaggio iniziale */
245 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)246 PiezoActuatorBeam2::InitialAssRes(SubVectorHandler& WorkVec,
247 const VectorHandler& XCurr)
248 {
249 return Beam2::InitialAssRes(WorkVec, XCurr);
250 }
251
252 /* PiezoActuatorBeam2 - end */
253
254
255 /* PiezoActuatorVEBeam2 - begin */
256
257 /* Funzioni di calcolo delle matrici */
258 void
AssStiffnessMat(FullSubMatrixHandler & WMA,FullSubMatrixHandler & WMB,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)259 PiezoActuatorVEBeam2::AssStiffnessMat(FullSubMatrixHandler& WMA,
260 FullSubMatrixHandler& WMB,
261 doublereal dCoef,
262 const VectorHandler& XCurr,
263 const VectorHandler& XPrimeCurr)
264 {
265 DEBUGCOUT("Entering PiezoActuatorVEBeam2::AssStiffnessMat()" << std::endl);
266
267 ViscoElasticBeam2::AssStiffnessMat(WMA, WMB, dCoef, XCurr, XPrimeCurr);
268
269 Mat3xN tmp1(iNumElec);
270 Mat3xN tmp2(iNumElec);
271
272 tmp1.LeftMult(R*dCoef, PiezoMat[STRAIN]);
273 WMA.Sub(1, 13, tmp1);
274 WMA.Add(7, 13, tmp1);
275
276 tmp2.LeftMult(Mat3x3(MatCross, p - pNode[NODE1]->GetXCurr()), tmp1);
277 WMA.Add(4, 13, tmp2);
278
279 tmp2.LeftMult(Mat3x3(MatCross, p - pNode[NODE2]->GetXCurr()), tmp1);
280 WMA.Sub(10, 13, tmp2);
281
282 tmp1.LeftMult(R*dCoef, PiezoMat[CURVAT]);
283 WMA.Sub(4, 13, tmp1);
284 WMA.Add(10, 13, tmp1);
285 }
286
287
288 void
AssStiffnessVec(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)289 PiezoActuatorVEBeam2::AssStiffnessVec(SubVectorHandler& WorkVec,
290 doublereal dCoef,
291 const VectorHandler& XCurr,
292 const VectorHandler& XPrimeCurr)
293 {
294 DEBUGCOUT("Entering PiezoActuatorVEBeam2::AssStiffnessVec()" << std::endl);
295
296 /*
297 * Riceve il vettore gia' dimensionato e con gli indici a posto
298 * per scrivere il residuo delle equazioni di equilibrio dei due nodi
299 */
300
301 /*
302 * Per la trattazione teorica, il riferimento e' il file ul-travi.tex
303 * (ora e' superato)
304 */
305
306 if (bFirstRes) {
307 /*
308 * AfterPredict ha gia' calcolato tutto;
309 * bFirstRes viene resettato da
310 * ViscoElasticBeam2::AssStiffnessVec
311 */
312
313 } else {
314 for (integer iCnt = 1; iCnt <= iNumElec; iCnt++) {
315 V.Put(iCnt, pvElecDofs[iCnt-1]->dGetX());
316 }
317 }
318
319 ViscoElasticBeam2::AssStiffnessVec(WorkVec, dCoef, XCurr, XPrimeCurr);
320 }
321
322
323 void
AddInternalForces(Vec6 & AzLoc)324 PiezoActuatorVEBeam2::AddInternalForces(Vec6& AzLoc)
325 {
326 AzLoc += Vec6(PiezoMat[STRAIN]*V, PiezoMat[CURVAT]*V);
327 }
328
329
330 /* Costruttore normale */
PiezoActuatorVEBeam2(unsigned int uL,const StructNode * pN1,const StructNode * pN2,const Vec3 & F1,const Vec3 & F2,const Mat3x3 & R1,const Mat3x3 & R2,const Mat3x3 & r,const ConstitutiveLaw6D * pD,int iEl,const ScalarDifferentialNode ** pEDof,const Mat3xN & Te,const Mat3xN & Tk,OrientationDescription ood,flag fOut)331 PiezoActuatorVEBeam2::PiezoActuatorVEBeam2(unsigned int uL,
332 const StructNode* pN1, const StructNode* pN2,
333 const Vec3& F1, const Vec3& F2,
334 const Mat3x3& R1, const Mat3x3& R2,
335 const Mat3x3& r,
336 const ConstitutiveLaw6D* pD,
337 int iEl,
338 const ScalarDifferentialNode **pEDof,
339 const Mat3xN& Te, const Mat3xN& Tk,
340 OrientationDescription ood,
341 flag fOut)
342 : Elem(uL, fOut),
343 ViscoElasticBeam2(uL, pN1, pN2, F1, F2, R1, R2, r, pD, ood, fOut),
344 iNumElec(iEl), pvElecDofs(pEDof), V(iEl)
345 {
346 ASSERT(iNumElec > 0);
347 ASSERT(pvElecDofs != NULL);
348
349 #ifdef DEBUG
350 for (int i = iNumElec; i-- > 0; ) {
351 ASSERT(pvElecDofs[i] != NULL);
352 }
353 #endif /* DEBUG */
354
355 PiezoMat[STRAIN].Copy(Te);
356 PiezoMat[CURVAT].Copy(Tk);
357 }
358
359
360 /* Distruttore banale */
~PiezoActuatorVEBeam2(void)361 PiezoActuatorVEBeam2::~PiezoActuatorVEBeam2(void)
362 {
363 ASSERT(pvElecDofs != NULL);
364 SAFEDELETEARR(pvElecDofs);
365 }
366
367
368 /* Contributo al file di restart */
369 std::ostream&
Restart(std::ostream & out) const370 PiezoActuatorVEBeam2::Restart(std::ostream& out) const
371 {
372 Restart_(out);
373 return out << "/* piezoelectric actuator NOT IMPLEMENTED YET */ ;"
374 << std::endl;
375 }
376
377
378 /*
379 * Dimensioni del workspace
380 */
381 void
WorkSpaceDim(integer * piNumRows,integer * piNumCols) const382 PiezoActuatorVEBeam2::WorkSpaceDim(integer* piNumRows,
383 integer* piNumCols) const
384 {
385 *piNumRows = 12;
386 *piNumCols = 12+iNumElec;
387 }
388
389
390 /* Settings iniziali, prima della prima soluzione */
391 void
SetValue(DataManager * pDM,VectorHandler & X,VectorHandler & XP,SimulationEntity::Hints * ph)392 PiezoActuatorVEBeam2::SetValue(DataManager *pDM,
393 VectorHandler& X, VectorHandler& XP,
394 SimulationEntity::Hints *ph)
395 {
396 /* se proprio non serve, si puo' eliminare */
397 ViscoElasticBeam2::SetValue(pDM, X, XP, ph);
398 }
399
400
401 /* Prepara i parametri di riferimento dopo la predizione */
402 void
AfterPredict(VectorHandler & X,VectorHandler & XP)403 PiezoActuatorVEBeam2::AfterPredict(VectorHandler& X, VectorHandler& XP)
404 {
405 /*
406 * Calcola le deformazioni, aggiorna il legame costitutivo
407 * e crea la FDE
408 */
409 for (integer iCnt = 1; iCnt <= iNumElec; iCnt++) {
410 V.Put(iCnt, pvElecDofs[iCnt-1]->dGetX());
411 }
412
413 ViscoElasticBeam2::AfterPredict(X, XP);
414 }
415
416
417 /* assemblaggio jacobiano */
418 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)419 PiezoActuatorVEBeam2::AssJac(VariableSubMatrixHandler& WorkMat,
420 doublereal dCoef,
421 const VectorHandler& XCurr,
422 const VectorHandler& XPrimeCurr)
423 {
424 DEBUGCOUT("Entering PiezoActuatorVEBeam2::AssJac();"
425 " will result in call to AssStiffnessMat()" << std::endl);
426
427 integer iNode1FirstMomIndex = pNode[NODE1]->iGetFirstMomentumIndex();
428 integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
429 integer iNode2FirstMomIndex = pNode[NODE2]->iGetFirstMomentumIndex();
430 integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
431
432 FullSubMatrixHandler& WM = WorkMat.SetFull();
433
434 /* Dimensiona la matrice, la azzera e pone gli indici corretti */
435 WM.ResizeReset(12, 12+iNumElec);
436
437 for (int iCnt = 1; iCnt <= 6; iCnt++) {
438 WM.PutRowIndex(iCnt, iNode1FirstMomIndex+iCnt);
439 WM.PutColIndex(iCnt, iNode1FirstPosIndex+iCnt);
440 WM.PutRowIndex(6+iCnt, iNode2FirstMomIndex+iCnt);
441 WM.PutColIndex(6+iCnt, iNode2FirstPosIndex+iCnt);
442 }
443
444 for (int iCnt = 1; iCnt <= iNumElec; iCnt++) {
445 WM.PutColIndex(12+iCnt,
446 pvElecDofs[iCnt-1]->iGetFirstColIndex()+1);
447 }
448
449 AssStiffnessMat(WM, WM, dCoef, XCurr, XPrimeCurr);
450
451 return WorkMat;
452 }
453
454
455 /* Contributo allo jacobiano durante l'assemblaggio iniziale */
456 VariableSubMatrixHandler&
InitialAssJac(VariableSubMatrixHandler & WorkMat,const VectorHandler & XCurr)457 PiezoActuatorVEBeam2::InitialAssJac(VariableSubMatrixHandler& WorkMat,
458 const VectorHandler& XCurr)
459 {
460 return ViscoElasticBeam2::InitialAssJac(WorkMat, XCurr);
461 }
462
463
464 /* Contributo al residuo durante l'assemblaggio iniziale */
465 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)466 PiezoActuatorVEBeam2::InitialAssRes(SubVectorHandler& WorkVec,
467 const VectorHandler& XCurr)
468 {
469 return ViscoElasticBeam2::InitialAssRes(WorkVec, XCurr);
470 }
471
472 /* PiezoActuatorVEBeam2 - end */
473
474