52 #include "Amesos_Klu.h" 56 #include "Epetra_BLAS.h" 57 #include "Epetra_CrsMatrix.h" 58 #include "Epetra_Export.h" 59 #include "Epetra_FECrsMatrix.h" 60 #include "Epetra_FEVector.h" 61 #include "Epetra_Import.h" 62 #include "Epetra_IntSerialDenseVector.h" 63 #include "Epetra_LAPACK.h" 64 #include "Epetra_Map.h" 65 #include "Epetra_MultiVector.h" 66 #include "Epetra_SerialDenseMatrix.h" 67 #include "Epetra_Vector.h" 68 #include "Teuchos_ParameterList.hpp" 69 #include "Teuchos_RCP.hpp" 70 #include "Teuchos_Assert.hpp" 71 #include "Teuchos_VerboseObject.hpp" 74 # include "Epetra_MpiComm.h" 77 # include "Epetra_SerialComm.h" 95 const double GLp_pi = 3.14159265358979323846;
105 int quadrature(
const int,
const int, Epetra_SerialDenseMatrix &,
106 Epetra_SerialDenseVector &);
110 Epetra_IntSerialDenseVector &,
111 Epetra_SerialDenseMatrix &,
112 Epetra_IntSerialDenseVector &,
113 Epetra_SerialDenseMatrix &,
114 Epetra_IntSerialDenseMatrix &,
115 Epetra_IntSerialDenseMatrix &,
116 const char geomFileBase[],
117 const bool trace =
true,
118 const bool dumpAll =
false 121 int lassembly(
const Epetra_SerialDenseMatrix &,
122 const Epetra_SerialDenseVector &,
123 const Epetra_SerialDenseMatrix &,
124 const Epetra_SerialDenseVector &,
125 const Epetra_SerialDenseVector &,
127 Epetra_SerialDenseMatrix &,
128 Epetra_SerialDenseVector &);
131 const Epetra_IntSerialDenseVector &,
132 const Epetra_SerialDenseMatrix &,
133 const Epetra_IntSerialDenseVector &,
134 const Epetra_SerialDenseMatrix &,
135 const Epetra_IntSerialDenseMatrix &,
136 const Epetra_IntSerialDenseMatrix &,
137 Teuchos::RCP<Epetra_FECrsMatrix> &,
138 Teuchos::RCP<Epetra_FEVector> &);
141 const Epetra_IntSerialDenseVector &,
142 const Epetra_SerialDenseMatrix &,
143 const Epetra_IntSerialDenseVector &,
144 const Epetra_SerialDenseMatrix &,
145 const Epetra_IntSerialDenseMatrix &,
146 const Epetra_IntSerialDenseMatrix &,
147 Teuchos::RCP<Epetra_FECrsMatrix> &,
148 Teuchos::RCP<Epetra_FEVector> &,
152 const Epetra_IntSerialDenseVector &,
153 const Epetra_SerialDenseMatrix &,
154 const Epetra_IntSerialDenseVector &,
155 const Epetra_SerialDenseMatrix &,
156 const Epetra_IntSerialDenseMatrix &,
157 const Epetra_IntSerialDenseMatrix &,
158 Teuchos::RCP<Epetra_FECrsMatrix> &,
159 Teuchos::RCP<Epetra_FECrsMatrix> &,
160 Teuchos::RCP<Epetra_FEVector> &);
163 const Epetra_Comm &Comm
164 ,
const Epetra_IntSerialDenseVector &ipindx
165 ,
const Epetra_SerialDenseMatrix &ipcoords
166 ,
const Epetra_IntSerialDenseVector &pindx
167 ,
const Epetra_SerialDenseMatrix &pcoords
168 ,
const Epetra_IntSerialDenseMatrix &t
169 ,
const Epetra_IntSerialDenseMatrix &e
170 ,Teuchos::RCP<Epetra_FECrsMatrix> *B
171 ,Teuchos::RCP<Epetra_FECrsMatrix> *R
175 const Epetra_IntSerialDenseVector &,
176 const Epetra_SerialDenseMatrix &,
177 const Epetra_IntSerialDenseVector &,
178 const Epetra_SerialDenseMatrix &,
179 const Epetra_IntSerialDenseMatrix &,
180 const Teuchos::RCP<const Epetra_MultiVector> &,
181 Teuchos::RCP<Epetra_FEVector> &);
184 const Epetra_IntSerialDenseVector &,
185 const Epetra_SerialDenseMatrix &,
186 const Epetra_IntSerialDenseVector &,
187 const Epetra_SerialDenseMatrix &,
188 const Epetra_IntSerialDenseMatrix &,
189 const Teuchos::RCP<const Epetra_MultiVector> &,
190 Teuchos::RCP<Epetra_FECrsMatrix> &);
193 const Epetra_IntSerialDenseVector &,
194 const Epetra_SerialDenseMatrix &,
195 const Epetra_IntSerialDenseVector &,
196 const Epetra_SerialDenseMatrix &,
197 const Epetra_IntSerialDenseMatrix &,
198 const Teuchos::RCP<const Epetra_MultiVector> &,
199 const Teuchos::RCP<const Epetra_MultiVector> &,
200 const Teuchos::RCP<const Epetra_MultiVector> &,
201 Teuchos::RCP<Epetra_FEVector> &);
203 int compproduct(Epetra_SerialDenseVector &,
double *,
double *);
205 int compproduct(Epetra_SerialDenseVector &,
double *,
double *,
double *);
207 double determinant(
const Epetra_SerialDenseMatrix &);
209 int inverse(
const Epetra_SerialDenseMatrix &, Epetra_SerialDenseMatrix &);
212 const int,
const int, Epetra_SerialDenseMatrix &,
213 Epetra_SerialDenseVector &);
215 void gpfctn(
const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv);
217 void g2pfctn(
const Epetra_SerialDenseVector & , Epetra_SerialDenseVector & );
219 void gfctn(
const Epetra_SerialDenseVector & , Epetra_SerialDenseVector & );
225 GLpYUEpetraDataPool::GLpYUEpetraDataPool(
226 Teuchos::RCP<const Epetra_Comm>
const& commptr
238 ipcoords_ = Teuchos::rcp(
new Epetra_SerialDenseMatrix() );
239 ipindx_ = Teuchos::rcp(
new Epetra_IntSerialDenseVector() );
240 pcoords_ = Teuchos::rcp(
new Epetra_SerialDenseMatrix() );
241 pindx_ = Teuchos::rcp(
new Epetra_IntSerialDenseVector() );
242 t_ = Teuchos::rcp(
new Epetra_IntSerialDenseMatrix() );
243 e_ = Teuchos::rcp(
new Epetra_IntSerialDenseMatrix() );
245 if( myfile && myfile[0] !=
'\0' ) {
246 meshreader(*commptr_,*ipindx_,*ipcoords_,*pindx_,*pcoords_,*t_,*e_,myfile,trace);
250 commptr_->NumProc(),commptr_->MyPID()
251 ,len_x,len_y,local_nx,local_ny,2
252 ,&*ipindx_,&*ipcoords_,&*pindx_,&*pcoords_,&*t_,&*e_
253 #ifdef GLPYUEPETRA_DATAPOOL_DUMP_ALL_MESH 254 ,&*Teuchos::VerboseObjectBase::getDefaultOStream(),
true 262 assemble( *commptr, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, *e_, A_, H_, b_ );
263 assemble_bdry( *commptr, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, *e_, &B_, &R_ );
266 Epetra_Map standardmap(A_->DomainMap());
267 q_ = Teuchos::rcp(
new Epetra_FEVector(standardmap,1));
268 int * qintvalues =
new int[standardmap.NumMyElements()];
269 double * qdvalues =
new double[standardmap.NumMyElements()];
270 standardmap.MyGlobalElements(qintvalues);
271 for (
int i = 0; i < standardmap.NumMyElements(); i++)
272 qdvalues[i]=cos( GLp_pi* ((*ipcoords_)(i,0)) ) * cos( GLp_pi* ((*ipcoords_)(i,1)) );
273 q_->ReplaceGlobalValues(standardmap.NumMyElements(), qintvalues, qdvalues);
274 q_->GlobalAssemble();
281 Teuchos::RCP<const Epetra_MultiVector> ey =
293 const Teuchos::RCP<const Epetra_MultiVector> & rhsy,
294 const Teuchos::RCP<const Epetra_MultiVector> & rhsu,
295 const Teuchos::RCP<const Epetra_MultiVector> & rhsp,
296 const Teuchos::RCP<Epetra_MultiVector> & y,
297 const Teuchos::RCP<Epetra_MultiVector> & u,
298 const Teuchos::RCP<Epetra_MultiVector> & p,
491 TEUCHOS_TEST_FOR_EXCEPT(
true);
532 Epetra_Map overlapmap(-1, pindx_->M(), (
int*)(pindx_)->A(), 1, *commptr_);
533 Epetra_Map standardmap(A_->DomainMap());
534 Teuchos::RCP<Epetra_MultiVector> yoverlap = Teuchos::rcp(
new Epetra_MultiVector(overlapmap, 1));
535 Epetra_Import Importer(overlapmap, standardmap);
536 yoverlap->Import(*y, Importer, Insert);
537 nonlinvec(*commptr_, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, yoverlap, Ny_);
543 Epetra_Map overlapmap(-1, pindx_->M(), (
int*)(pindx_)->A(), 1, *commptr_);
544 Epetra_Map standardmap(A_->DomainMap());
545 Teuchos::RCP<Epetra_MultiVector> yoverlap = Teuchos::rcp(
new Epetra_MultiVector(overlapmap, 1));
546 Epetra_Import Importer(overlapmap, standardmap);
547 yoverlap->Import(*y, Importer, Insert);
548 nonlinjac(*commptr_, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, yoverlap, Npy_);
554 Epetra_Map standardmap(A_->DomainMap());
555 Epetra_Map bdryctrlmap(B_->DomainMap());
559 int numstates = standardmap.NumGlobalElements();
561 int nummystates = standardmap.NumMyElements();
562 int nummycontrols = bdryctrlmap.NumMyElements();
564 Epetra_IntSerialDenseVector KKTmapindx(2*nummystates+nummycontrols);
567 Epetra_IntSerialDenseVector states(nummystates);
568 Epetra_IntSerialDenseVector controls(nummycontrols);
569 standardmap.MyGlobalElements(states.Values());
570 bdryctrlmap.MyGlobalElements(controls.Values());
571 for (
int i=0; i<nummystates; i++) {
572 KKTmapindx(i) = states(i);
573 KKTmapindx(nummystates+nummycontrols+i) = 2*numstates+states(i);
575 for (
int i=0; i<nummycontrols; i++) {
576 KKTmapindx(nummystates+i) = numstates+controls(i);
578 Epetra_Map KKTmap(-1, 2*nummystates+nummycontrols, KKTmapindx.Values(), indexBase, *(commptr_));
583 Augmat_ = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, KKTmap, 0));
587 for (
int i=0; i<nummystates+nummycontrols; i++) {
588 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], 1, one, KKTmapindx.Values()+i);
593 Epetra_SerialDenseVector values(nummyentries);
594 Epetra_IntSerialDenseVector indices(nummyentries);
596 for (
int i=0; i<nummystates; i++) {
597 nummyentries = A_->NumMyEntries(i);
598 values.Resize(nummyentries);
599 indices.Resize(nummyentries);
600 A_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
602 if (nummyentries > 0)
603 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(),
605 nummyentries = Npy_->NumMyEntries(i);
606 values.Resize(nummyentries);
607 indices.Resize(nummyentries);
608 Npy_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
610 if (nummyentries > 0)
611 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(),
615 for (
int i=0; i<nummystates; i++) {
616 nummyentries = B_->NumMyEntries(i);
617 values.Resize(nummyentries);
618 indices.Resize(nummyentries);
619 B_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
621 for (
int j=0; j<nummyentries; j++)
622 indices[j] = indices[j]+numstates;
623 if (nummyentries > 0)
624 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(),
629 Epetra_CrsMatrix & transA =
dynamic_cast<Epetra_CrsMatrix&
>(transposer(*A_));
630 Epetra_CrsMatrix & transB =
dynamic_cast<Epetra_CrsMatrix&
>(transposer(*B_));
631 Epetra_CrsMatrix & transNpy =
dynamic_cast<Epetra_CrsMatrix&
>(transposer(*Npy_));
633 for (
int i=0; i<nummystates; i++) {
634 nummyentries = transA.NumMyEntries(i);
635 values.Resize(nummyentries);
636 indices.Resize(nummyentries);
637 transA.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
639 for (
int j=0; j<nummyentries; j++)
640 indices[j] = indices[j]+2*numstates;
641 if (nummyentries > 0)
642 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], nummyentries, values.Values(),
644 nummyentries = transNpy.NumMyEntries(i);
645 values.Resize(nummyentries);
646 indices.Resize(nummyentries);
647 transNpy.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
649 for (
int j=0; j<nummyentries; j++)
650 indices[j] = indices[j]+2*numstates;
651 if (nummyentries > 0)
652 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], nummyentries, values.Values(),
656 for (
int i=0; i<nummystates; i++) {
657 nummyentries = transB.NumMyEntries(i);
658 values.Resize(nummyentries);
659 indices.Resize(nummyentries);
660 transB.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
662 for (
int j=0; j<nummyentries; j++)
663 indices[j] = indices[j]+2*numstates;
665 if (nummyentries > 0)
666 err = Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+numstates, nummyentries,
667 values.Values(), indices.Values());
670 std::cout <<
"Insertion of entries failed:\n";
671 std::cout << indices;
672 std::cout << nummyentries << std::endl;
673 std::cout <<
"at row: " << KKTmapindx.Values()[i]+numstates << std::endl << std::endl;
677 Augmat_->FillComplete(KKTmap, KKTmap);
689 Epetra_SerialDenseVector product(4);
696 N.Shape(Nodes.M(),3);
697 for (
int i=0; i<Nodes.M(); i++) {
698 N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
703 Nx1.Shape(Nodes.M(),3);
704 for (
int i=0; i<Nodes.M(); i++) {
710 Nx2.Shape(Nodes.M(),3);
711 for (
int i=0; i<Nodes.M(); i++) {
720 S1(0,0)= 1.0; S1(0,1)=-1.0; S1(0,2)= 0.0;
721 S1(1,0)=-1.0; S1(1,1)= 1.0; S1(1,2)= 0.0;
722 S1(2,0)= 0.0; S1(2,1)= 0.0; S1(2,2)= 0.0;
724 S2(0,0)= 2.0; S2(0,1)=-1.0; S2(0,2)=-1.0;
725 S2(1,0)=-1.0; S2(1,1)= 0.0; S2(1,2)= 1.0;
726 S2(2,0)=-1.0; S2(2,1)= 1.0; S2(2,2)= 0.0;
728 S3(0,0)= 1.0; S3(0,1)= 0.0; S3(0,2)=-1.0;
729 S3(1,0)= 0.0; S3(1,1)= 0.0; S3(1,2)= 0.0;
730 S3(2,0)=-1.0; S3(2,1)= 0.0; S3(2,2)= 1.0;
734 Nw.Multiply(
'T',
'N', 1.0, N, Weights, 0.0);
738 for (
int i=0; i<3; i++) {
739 for (
int j=0; j<3; j++) {
741 NNw(i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
746 NNNw =
new Epetra_SerialDenseMatrix[3];
747 NNNw[0].Shape(3,3); NNNw[1].Shape(3,3); NNNw[2].Shape(3,3);
748 for (
int i=0; i<3; i++) {
749 for (
int j=0; j<3; j++) {
750 for (
int k=0; k<3; k++) {
752 NNNw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
758 NdNdx1Nw =
new Epetra_SerialDenseMatrix[3];
759 NdNdx1Nw[0].Shape(3,3); NdNdx1Nw[1].Shape(3,3); NdNdx1Nw[2].Shape(3,3);
760 for (
int i=0; i<3; i++) {
761 for (
int j=0; j<3; j++) {
762 for (
int k=0; k<3; k++) {
764 NdNdx1Nw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
770 NdNdx2Nw =
new Epetra_SerialDenseMatrix[3];
771 NdNdx2Nw[0].Shape(3,3); NdNdx2Nw[1].Shape(3,3); NdNdx2Nw[2].Shape(3,3);
772 for (
int i=0; i<3; i++) {
773 for (
int j=0; j<3; j++) {
774 for (
int k=0; k<3; k++) {
776 NdNdx2Nw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
785 os <<
"\n\nQuadrature nodes:";
786 os <<
"\n-----------------";
788 os <<
"\n\nQuadrature weights:";
789 os <<
"\n-------------------\n";
791 os <<
"\n\nNodal basis functions:";
792 os <<
"\n----------------------";
794 os <<
"\n\nIntegrals of combinations of partial derivatives:";
795 os <<
"\n-------------------------------------------------";
799 os <<
"\n\nIntegrals of basis functions:";
800 os <<
"\n-----------------------------\n";
802 os <<
"\n\nIntegrals of products N(i)*N(j):";
803 os <<
"\n--------------------------------\n";
805 os <<
"\n\nIntegrals of products N(i)*N(j)*N(k):";
806 os <<
"\n-------------------------------------\n";
807 NNNw[0].Print(os); NNNw[1].Print(os); NNNw[2].Print(os);
808 os <<
"\n\nIntegrals of products N(i)*(dN(j)/dx1)*N(k):";
809 os <<
"\n--------------------------------------------\n";
810 NdNdx1Nw[0].Print(os); NdNdx1Nw[1].Print(os); NdNdx1Nw[2].Print(os);
811 os <<
"\n\nIntegrals of products N(i)*(dN(j)/dx2)*N(k):";
812 os <<
"\n--------------------------------------------\n";
813 NdNdx2Nw[0].Print(os); NdNdx2Nw[1].Print(os); NdNdx2Nw[2].Print(os);
829 double *first,
double *second)
831 for (
int i=0; i<product.M(); i++) {
832 product[i] = first[i]*second[i];
838 double *first,
double *second,
double *third)
840 for (
int i=0; i<product.M(); i++) {
841 product[i] = first[i]*second[i]*third[i];
891 const Epetra_Comm &Comm
892 ,
const Epetra_IntSerialDenseVector &ipindx
893 ,
const Epetra_SerialDenseMatrix &ipcoords
894 ,
const Epetra_IntSerialDenseVector &pindx
895 ,
const Epetra_SerialDenseMatrix &pcoords
896 ,
const Epetra_IntSerialDenseMatrix &t
897 ,
const Epetra_IntSerialDenseMatrix &e
898 ,Teuchos::RCP<Epetra_FECrsMatrix> *B_out
899 ,Teuchos::RCP<Epetra_FECrsMatrix> *R_out
905 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY 906 Teuchos::RCP<Teuchos::FancyOStream>
907 out = Teuchos::VerboseObjectBase::getDefaultOStream();
908 Teuchos::OSTab tab(out);
909 *out <<
"\nEntering assemble_bdry(...) ...\n";
912 int numLocElems = t.M();
913 int numLocEdges = e.M();
921 Epetra_IntSerialDenseVector BdryNode(2*numLocEdges);
922 for (
int j=0; j<numLocEdges; j++) {
923 BdryNode[j] = e(j,0);
924 BdryNode[j+numLocEdges] = e(j,1);
926 std::sort(BdryNode.Values(), BdryNode.Values()+2*numLocEdges);
927 lastindx = std::unique(BdryNode.Values(), BdryNode.Values()+2*numLocEdges);
928 const int numBdryNodes = lastindx - BdryNode.Values();
929 BdryNode.Resize(numBdryNodes);
934 Epetra_IntSerialDenseVector MyBdryNode(numBdryNodes);
935 lastindx = std::set_intersection(
936 BdryNode.Values(), BdryNode.Values()+numBdryNodes,
937 ipindx.Values(), ipindx.Values()+ipindx.M(),
940 const int numMyBdryNodes = lastindx - MyBdryNode.Values();
941 MyBdryNode.Resize(numMyBdryNodes);
946 Epetra_Map standardmap(-1, ipindx.M(),
const_cast<int*
>(ipindx.A()), indexBase, Comm);
947 Epetra_Map overlapmap(-1, pindx.M(),
const_cast<int*
>(pindx.A()), indexBase, Comm);
948 Epetra_Map mybdryctrlmap(-1, numMyBdryNodes, const_cast<int*>(MyBdryNode.A()), indexBase, Comm);
953 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY 954 *out <<
"\nstandardmap:\n";
955 standardmap.Print(Teuchos::OSTab(out).o());
956 *out <<
"\nmybdryctrlmap:\n";
957 mybdryctrlmap.Print(Teuchos::OSTab(out).o());
963 Teuchos::RCP<Epetra_FECrsMatrix>
964 B = rcp(
new Epetra_FECrsMatrix(Copy,standardmap,0)),
965 R = rcp(
new Epetra_FECrsMatrix(Copy,standardmap,0));
969 const int numNodesPerEdge = 2;
970 Epetra_IntSerialDenseVector epetra_nodes(numNodesPerEdge);
976 Epetra_SerialDenseMatrix Bt(2,2);
979 for(
int i=0; i < numLocEdges; i++ ) {
982 global_id_0 = e(i,0),
983 global_id_1 = e(i,1),
984 local_id_0 = overlapmap.LID(global_id_0),
985 local_id_1 = overlapmap.LID(global_id_1);
987 epetra_nodes(0) = global_id_0;
988 epetra_nodes(1) = global_id_1;
991 x0 = pcoords(local_id_0,0),
992 y0 = pcoords(local_id_0,1),
993 x1 = pcoords(local_id_1,0),
994 y1 = pcoords(local_id_1,1);
996 const double l = sqrt(pow(x0-x1,2) + pow(y0-y1,2));
999 const double l_sixth = l * (1.0/6.0);
1000 Bt(0,0) = l_sixth * 2.0;
1001 Bt(0,1) = l_sixth * 1.0;
1002 Bt(1,0) = l_sixth * 1.0;
1003 Bt(1,1) = l_sixth * 2.0;
1005 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY 1007 <<
"\nEdge global nodes = ("<<global_id_0<<
","<<global_id_1<<
")," 1008 <<
" local nodes = ("<<local_id_0<<
","<<local_id_1<<
")," 1009 <<
" (x0,y0)(x1,y1)=("<<x0<<
","<<y0<<
")("<<x1<<
","<<y1<<
")," 1010 <<
" Bt = ["<<Bt(0,0)<<
","<<Bt(0,1)<<
";"<<Bt(1,0)<<
","<<Bt(1,1)<<
"]\n";
1013 const int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
1014 err = B->InsertGlobalValues(epetra_nodes,Bt,format);
1015 if (err<0)
return(err);
1016 err = R->InsertGlobalValues(epetra_nodes,Bt,format);
1017 if (err<0)
return(err);
1035 err = B->GlobalAssemble(mybdryctrlmap,standardmap,
true);
1036 if (err<0)
return(err);
1037 err = R->GlobalAssemble(mybdryctrlmap,mybdryctrlmap,
true);
1038 if (err<0)
return(err);
1040 if(B_out) *B_out = B;
1041 if(R_out) *R_out = R;
1043 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY 1045 B->Print(Teuchos::OSTab(out).o());
1046 *out <<
"\nLeaving assemble_bdry(...) ...\n";
1110 const Epetra_IntSerialDenseVector & ipindx,
1111 const Epetra_SerialDenseMatrix & ipcoords,
1112 const Epetra_IntSerialDenseVector & pindx,
1113 const Epetra_SerialDenseMatrix & pcoords,
1114 const Epetra_IntSerialDenseMatrix & t,
1115 const Epetra_IntSerialDenseMatrix & e,
1116 RCP<Epetra_FECrsMatrix> & A,
1117 RCP<Epetra_FECrsMatrix> & M,
1118 RCP<Epetra_FEVector> & b)
1121 int myPID = Comm.MyPID();
1122 int numProcs = Comm.NumProc();
1125 int numLocElems = t.M();
1126 int numNodesPerElem = 3;
1130 Epetra_Map standardmap(-1, ipindx.M(), (
int*)ipindx.A(), indexBase, Comm);
1131 Epetra_Map overlapmap(-1, pindx.M(), (
int*)pindx.A(), indexBase, Comm);
1133 int* nodes =
new int[numNodesPerElem];
1134 int i=0, j=0, err=0;
1136 A = rcp(
new Epetra_FECrsMatrix(Copy, standardmap, 0));
1137 M = rcp(
new Epetra_FECrsMatrix(Copy, standardmap, 0));
1138 b = rcp(
new Epetra_FEVector(standardmap,1));
1141 int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
1142 Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
1146 Epetra_SerialDenseMatrix At;
1147 Epetra_SerialDenseVector bt;
1148 Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
1150 Epetra_SerialDenseVector k(numNodesPerElem);
1151 for (i=0; i< numNodesPerElem; i++) k(i)=1.0;
1152 Epetra_SerialDenseMatrix c(numNodesPerElem,2);
1153 for (i=0; i< numNodesPerElem; i++) { c(i,0)=0.0; c(i,1)=0.0; }
1154 Epetra_SerialDenseVector r(numNodesPerElem);
1155 for (i=0; i< numNodesPerElem; i++) r(i)=0.0;
1156 Epetra_SerialDenseVector f(numNodesPerElem);
1157 for (i=0; i< numNodesPerElem; i++) f(i)=0.0;
1158 Epetra_SerialDenseVector g(2);
1160 Epetra_SerialDenseVector sigma(2);
1161 sigma(0)=0.0; sigma(1)=0.0;
1162 for(i=0; i<numLocElems; i++) {
1163 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
1164 for (j=0; j<numNodesPerElem; j++) {
1166 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
1167 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
1169 f(j) = cos(
GLp_pi*vertices(j,0))*cos(
GLp_pi*vertices(j,1)) *
1172 lassembly(vertices, k, c, r, f, usr_par, At, bt);
1173 err = A->InsertGlobalValues(epetra_nodes, At, format);
1174 if (err<0)
return(err);
1175 err = b->SumIntoGlobalValues(epetra_nodes, bt);
1176 if (err<0)
return(err);
1182 Epetra_IntSerialDenseMatrix neumann(e.M(),2);
1184 for (i=0; i<e.M(); i++) {
1186 neumann(j,0) = e(i,0); neumann(j,1) = e(i,1);
1190 neumann.Reshape(j,2);
1192 if (neumann.M() != 0) {
1194 Epetra_SerialDenseMatrix quadnodes;
1195 Epetra_SerialDenseVector quadweights;
1196 Epetra_SerialDenseMatrix N;
1197 Epetra_SerialDenseMatrix NN;
1198 Epetra_SerialDenseVector product(2);
1205 N.Shape(quadnodes.M(),2);
1206 for (i=0; i<quadnodes.M(); i++) {
1207 N(i,0) = 1.0 - quadnodes(i,0);
1208 N(i,1) = quadnodes(i,0);
1213 for (i=0; i<2; i++) {
1214 for (j=0; j<2; j++) {
1216 NN(i,j) = blas.DOT(quadweights.M(), quadweights.A(), product.A());
1220 Epetra_IntSerialDenseVector neumann_nodes(2);
1221 Epetra_SerialDenseVector neumann_add(2);
1223 for (i=0; i<neumann.M(); i++) {
1224 neumann_nodes(0) = neumann(i,0); neumann_nodes(1) = neumann(i,1);
1225 neumann_add(0) = pcoords(overlapmap.LID(neumann_nodes(0)),0)
1226 -pcoords(overlapmap.LID(neumann_nodes(1)),0);
1227 neumann_add(1) = pcoords(overlapmap.LID(neumann_nodes(0)),1)
1228 -pcoords(overlapmap.LID(neumann_nodes(1)),1);
1229 l = blas.NRM2(neumann_add.M(), neumann_add.A());
1230 neumann_add.Multiply(
'N',
'N', 1.0, NN, g, 0.0);
1231 neumann_add.Scale(l);
1232 err = b->SumIntoGlobalValues(neumann_nodes, neumann_add);
1233 if (err<0)
return(err);
1238 Epetra_IntSerialDenseMatrix robin(e.M(),2);
1240 for (i=0; i<e.M(); i++) {
1242 robin(j,0) = e(i,0); robin(j,1) = e(i,1);
1248 if (robin.M() != 0) {
1250 Epetra_SerialDenseMatrix quadnodes;
1251 Epetra_SerialDenseVector quadweights;
1252 Epetra_SerialDenseMatrix N;
1253 Epetra_SerialDenseMatrix NN;
1254 Epetra_SerialDenseMatrix * NNN;
1255 Epetra_SerialDenseVector product(2);
1262 N.Shape(quadnodes.M(),2);
1263 for (i=0; i<quadnodes.M(); i++) {
1264 N(i,0) = 1.0 - quadnodes(i,0);
1265 N(i,1) = quadnodes(i,0);
1270 for (i=0; i<2; i++) {
1271 for (j=0; j<2; j++) {
1273 NN(i,j) = blas.DOT(quadweights.M(), quadweights.A(), product.A());
1278 NNN =
new Epetra_SerialDenseMatrix[2];
1279 NNN[0].Shape(2,2); NNN[1].Shape(2,2);
1280 for (i=0; i<2; i++) {
1281 for (j=0; j<2; j++) {
1282 for (
int k=0; k<2; k++) {
1284 NNN[k](i,j) = blas.DOT(quadweights.M(), quadweights.A(),
1290 Epetra_IntSerialDenseVector robin_nodes(2);
1291 Epetra_SerialDenseVector robin_b_add(2);
1292 Epetra_SerialDenseMatrix robin_A_add(2,2);
1294 for (i=0; i<robin.M(); i++) {
1295 robin_nodes(0) = robin(i,0); robin_nodes(1) = robin(i,1);
1297 robin_b_add(0) = pcoords(overlapmap.LID(robin_nodes(0)),0)
1298 -pcoords(overlapmap.LID(robin_nodes(1)),0);
1299 robin_b_add(1) = pcoords(overlapmap.LID(robin_nodes(0)),1)
1300 -pcoords(overlapmap.LID(robin_nodes(1)),1);
1301 l = blas.NRM2(robin_b_add.M(), robin_b_add.A());
1302 robin_b_add.Multiply(
'N',
'N', 1.0, NN, g, 0.0);
1303 robin_b_add.Scale(l);
1304 err = b->SumIntoGlobalValues(robin_nodes, robin_b_add);
1305 if (err<0)
return(err);
1307 NNN[0].Scale(sigma(0)); NNN[1].Scale(sigma(1));
1308 robin_A_add += NNN[0]; robin_A_add += NNN[1];
1309 robin_A_add.Scale(l);
1310 err = A->InsertGlobalValues(robin_nodes, robin_A_add, format);
1311 if (err<0)
return(err);
1318 Epetra_IntSerialDenseMatrix dirichlet(e.M(),2);
1320 for (i=0; i<e.M(); i++) {
1322 dirichlet(j,0) = e(i,0); dirichlet(j,1) = e(i,1);
1326 dirichlet.Reshape(j,2);
1335 Epetra_SerialDenseMatrix Mt;
1337 for (i=0; i< numNodesPerElem; i++) k(i)=0.0;
1338 for (i=0; i< numNodesPerElem; i++) { c(i,0)=0.0; c(i,1)=0.0; }
1339 for (i=0; i< numNodesPerElem; i++) r(i)=1.0;
1340 for (i=0; i< numNodesPerElem; i++) f(i)=0.0;
1342 sigma(0)=0.0; sigma(1)=0.0;
1343 for(i=0; i<numLocElems; i++) {
1344 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
1345 for (j=0; j<numNodesPerElem; j++) {
1346 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
1347 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
1349 lassembly(vertices, k, c, r, f, usr_par, Mt, bt);
1350 err = M->InsertGlobalValues(epetra_nodes, Mt, format);
1351 if (err<0)
return(err);
1360 err = A->GlobalAssemble();
1361 if (err<0)
return(err);
1363 err = b->GlobalAssemble();
1364 if (err<0)
return(err);
1366 err = M->GlobalAssemble();
1367 if (err<0)
return(err);
1415 const Epetra_SerialDenseVector & k,
1416 const Epetra_SerialDenseMatrix & c,
1417 const Epetra_SerialDenseVector & r,
1418 const Epetra_SerialDenseVector & f,
1420 Epetra_SerialDenseMatrix & At,
1421 Epetra_SerialDenseVector & bt)
1424 Epetra_SerialDenseMatrix B(2,2);
1425 Epetra_SerialDenseMatrix b(2,2);
1426 Epetra_SerialDenseMatrix BtB(2,2);
1427 Epetra_SerialDenseMatrix C(2,2);
1428 Epetra_SerialDenseMatrix M1(3,3);
1429 Epetra_SerialDenseMatrix M2(3,3);
1430 Epetra_SerialDenseMatrix M3(3,3);
1431 Epetra_SerialDenseMatrix tmp(3,3);
1437 for(
int i=0; i<2; i++) {
1438 B(i,0) = vertices(1,i)-vertices(0,i);
1439 B(i,1) = vertices(2,i)-vertices(0,i);
1445 BtB.Multiply(
'T',
'N', 1.0, B, B, 0.0);
1451 tmp = usr_par.
S1; tmp.Scale(C(0,0));
1453 tmp = usr_par.
S2; tmp.Scale(C(0,1));
1455 tmp = usr_par.
S3; tmp.Scale(C(1,1));
1457 M1.Scale( (k(0)*usr_par.
Nw(0) + k(1)*usr_par.
Nw(1) +
1458 k(2)*usr_par.
Nw(2)) * adB );
1461 tmp = usr_par.
NdNdx1Nw[0]; tmp.Scale(b(0,0)*c(0,0)+b(0,1)*c(0,1));
1463 tmp = usr_par.
NdNdx1Nw[1]; tmp.Scale(b(0,0)*c(1,0)+b(0,1)*c(1,1));
1465 tmp = usr_par.
NdNdx1Nw[2]; tmp.Scale(b(0,0)*c(2,0)+b(0,1)*c(2,1));
1467 tmp = usr_par.
NdNdx2Nw[0]; tmp.Scale(b(1,0)*c(0,0)+b(1,1)*c(0,1));
1469 tmp = usr_par.
NdNdx2Nw[1]; tmp.Scale(b(1,0)*c(1,0)+b(1,1)*c(1,1));
1471 tmp = usr_par.
NdNdx2Nw[2]; tmp.Scale(b(1,0)*c(2,0)+b(1,1)*c(2,1));
1476 tmp = usr_par.
NNNw[0]; tmp.Scale(r(0));
1478 tmp = usr_par.
NNNw[1]; tmp.Scale(r(1));
1480 tmp = usr_par.
NNNw[2]; tmp.Scale(r(2));
1488 bt.Multiply(
'T',
'N', adB, usr_par.
NNw, f, 0.0);
1501 Epetra_SerialDenseMatrix & inv)
1503 Epetra_LAPACK lapack;
1506 Epetra_IntSerialDenseVector ipiv(dim);
1507 Epetra_SerialDenseVector work(2*dim);
1512 lapack.GETRF(dim, dim, inv.A(), dim, ipiv.A(), &info);
1513 lapack.GETRI(dim, inv.A(), dim, ipiv.A(), work.A(), &dim, &info);
1527 Epetra_LAPACK lapack;
1532 Epetra_IntSerialDenseVector ipiv(dim);
1534 Epetra_SerialDenseMatrix mymat(mat);
1535 lapack.GETRF(dim, dim, mymat.A(), dim, ipiv.A(), &info);
1538 for (
int i=0; i<dim; i++) {
1543 for (
int i=0; i<dim; i++) {
1544 if ((ipiv[i]-1) != i)
1548 det *= pow((
double)(-1.0),swaps);
1554 Epetra_IntSerialDenseVector & ipindx,
1555 Epetra_SerialDenseMatrix & ipcoords,
1556 Epetra_IntSerialDenseVector & pindx,
1557 Epetra_SerialDenseMatrix & pcoords,
1558 Epetra_IntSerialDenseMatrix & t,
1559 Epetra_IntSerialDenseMatrix & e,
1560 const char geomFileBase[],
1565 int MyPID = Comm.MyPID();
1567 const int FileNameSize = 120;
1568 char FileName[FileNameSize];
1569 TEUCHOS_TEST_FOR_EXCEPT(static_cast<int>(std::strlen(geomFileBase) + 5) > FileNameSize);
1570 sprintf(FileName,
"%s.%03d", geomFileBase, MyPID);
1573 std::ifstream file_in(FileName);
1574 TEUCHOS_TEST_FOR_EXCEPTION(
1575 file_in.eof(), std::logic_error
1576 ,
"Error, the file \""<<FileName<<
"\" could not be opened for input!" 1581 fid = fopen(FileName,
"r");
1583 if(trace) printf(
"\nReading node info from %s ...\n", FileName);
1584 int numip = 0, numcp = 0;
1585 fscanf(fid,
"%d %d", &numip, &numcp);
1586 const int nump = numip + numcp;
1587 if(trace) printf(
"\nnumip = %d, numcp = %d, nump = %d\n", numip, numcp, nump );
1589 ipcoords.Shape(numip, 2);
1591 pcoords.Shape(nump, 2);
1592 for (
int i=0; i<numip; i++) {
1593 fscanf(fid,
"%d %lf %lf %*d",&ipindx(i),&ipcoords(i,0),&ipcoords(i,1));
1594 if(trace&&dumpAll) printf(
"%d %lf %lf\n",ipindx(i),ipcoords(i,0),ipcoords(i,1));
1595 pindx(i) = ipindx(i);
1596 pcoords(i,0) = ipcoords(i,0); pcoords(i,1) = ipcoords(i,1);
1598 for (
int i=numip; i<nump; i++) {
1599 fscanf(fid,
"%d %lf %lf %*d",&pindx(i),&pcoords(i,0),&pcoords(i,1));
1602 fscanf(fid,
"%*[^\n]");
1603 fscanf(fid,
"%*1[\n]");
1605 fscanf(fid,
"%*[^\n]");
1606 fscanf(fid,
"%*1[\n]");
1608 for (
int i=0; i<nump; i++) {
1609 fscanf(fid,
"%*[^\n]");
1610 fscanf(fid,
"%*1[\n]");
1613 if(trace) printf(
"\nReading element info from %s ...\n", FileName);
1615 fscanf(fid,
"%d", &numelems);
1616 if(trace) printf(
"\nnumelems = %d\n", numelems );
1617 t.Shape(numelems,3);
1618 for (
int i=0; i<numelems; i++) {
1619 fscanf(fid,
"%d %d %d", &t(i,0), &t(i,1), &t(i,2));
1620 if(trace&&dumpAll) printf(
"%d %d %d\n", t(i,0), t(i,1), t(i,2));
1623 if(trace) printf(
"\nReading boundry edge info from %s ...\n", FileName);
1625 fscanf(fid,
"%d",&numedges);
1626 if(trace) printf(
"\nnumedges = %d\n", numedges );
1627 e.Shape(numedges,3);
1628 for (
int i=0; i<numedges; i++) {
1629 fscanf(fid,
"%d %d %d", &e(i,0), &e(i,1), &e(i,2));
1630 if(trace&&dumpAll) printf(
"%d %d %d\n", e(i,0), e(i,1), e(i,2));
1634 if(trace) printf(
"Done reading mesh.\n");
1677 const Epetra_IntSerialDenseVector & ipindx,
1678 const Epetra_SerialDenseMatrix & ipcoords,
1679 const Epetra_IntSerialDenseVector & pindx,
1680 const Epetra_SerialDenseMatrix & pcoords,
1681 const Epetra_IntSerialDenseMatrix & t,
1682 const Teuchos::RCP<const Epetra_MultiVector> & y,
1683 const Teuchos::RCP<const Epetra_MultiVector> & s,
1684 const Teuchos::RCP<const Epetra_MultiVector> & lambda,
1685 Teuchos::RCP<Epetra_FEVector> & hessvec)
1688 int myPID = Comm.MyPID();
1689 int numProcs = Comm.NumProc();
1691 int numLocNodes = pindx.M();
1692 int numMyLocNodes = ipindx.M();
1693 int numLocElems = t.M();
1694 int numNodesPerElem = 3;
1698 Epetra_Map standardmap(-1, numMyLocNodes, (
int*)ipindx.A(), indexBase, Comm);
1699 Epetra_Map overlapmap(-1, numLocNodes, (
int*)pindx.A(), indexBase, Comm);
1701 hessvec = rcp(
new Epetra_FEVector(standardmap,1));
1703 int* nodes =
new int[numNodesPerElem];
1704 int i=0, j=0, err=0;
1707 Epetra_SerialDenseMatrix Nodes;
1708 Epetra_SerialDenseVector Weights;
1710 int numQuadPts = Nodes.M();
1714 Epetra_SerialDenseMatrix N;
1715 N.Shape(numQuadPts,3);
1716 for (
int i=0; i<numQuadPts; i++) {
1717 N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
1718 N(i,1) = Nodes(i,0);
1719 N(i,2) = Nodes(i,1);
1723 Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
1724 Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
1726 Epetra_SerialDenseVector ly;
1727 Epetra_SerialDenseVector Nly;
1728 Epetra_SerialDenseVector ls;
1729 Epetra_SerialDenseVector Nls;
1730 Epetra_SerialDenseVector llambda;
1731 Epetra_SerialDenseVector Nllambda;
1732 Epetra_SerialDenseVector lgfctn;
1733 Epetra_SerialDenseVector lgfctnNi;
1734 Epetra_SerialDenseVector lgfctnNls;
1735 Epetra_SerialDenseVector lhessvec;
1737 ly.Size(numNodesPerElem);
1738 Nly.Size(numQuadPts);
1739 ls.Size(numNodesPerElem);
1740 Nls.Size(numQuadPts);
1741 llambda.Size(numNodesPerElem);
1742 Nllambda.Size(numQuadPts);
1743 lgfctn.Size(numQuadPts);
1744 lgfctnNi.Size(numQuadPts);
1745 lgfctnNls.Size(numQuadPts);
1746 lhessvec.Size(numNodesPerElem);
1748 Epetra_SerialDenseMatrix B(2,2);
1751 for(i=0; i<numLocElems; i++) {
1753 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
1754 for (j=0; j<numNodesPerElem; j++) {
1755 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
1756 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
1760 for(
int i=0; i<2; i++) {
1761 B(i,0) = vertices(1,i)-vertices(0,i);
1762 B(i,1) = vertices(2,i)-vertices(0,i);
1767 for (j=0; j<numNodesPerElem; j++) {
1768 ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
1769 ls(j) = (*((*s)(0)))[overlapmap.LID(nodes[j])];
1770 llambda(j) = (*((*lambda)(0)))[overlapmap.LID(nodes[j])];
1773 Nly.Multiply(
'N',
'N', 1.0, N, ly, 0.0);
1774 Nls.Multiply(
'N',
'N', 1.0, N, ls, 0.0);
1775 Nllambda.Multiply(
'N',
'N', 1.0, N, llambda, 0.0);
1778 for (
int i=0; i<numNodesPerElem; i++) {
1780 compproduct(lgfctnNls, lgfctnNi.A(), Nls.A(), Nllambda.A());
1781 lhessvec(i) = adB*lgfctnNls.Dot(Weights);
1784 err = hessvec->SumIntoGlobalValues(epetra_nodes, lhessvec);
1785 if (err<0)
return(err);
1790 err = hessvec->GlobalAssemble();
1791 if (err<0)
return(err);
1803 for (
int i=0; i<v.M(); i++) {
1841 const Epetra_IntSerialDenseVector & ipindx,
1842 const Epetra_SerialDenseMatrix & ipcoords,
1843 const Epetra_IntSerialDenseVector & pindx,
1844 const Epetra_SerialDenseMatrix & pcoords,
1845 const Epetra_IntSerialDenseMatrix & t,
1846 const Teuchos::RCP<const Epetra_MultiVector> & y,
1847 Teuchos::RCP<Epetra_FECrsMatrix> & Gp)
1850 int myPID = Comm.MyPID();
1851 int numProcs = Comm.NumProc();
1853 int numLocNodes = pindx.M();
1854 int numMyLocNodes = ipindx.M();
1855 int numLocElems = t.M();
1856 int numNodesPerElem = 3;
1860 Epetra_Map standardmap(-1, numMyLocNodes, (
int*)ipindx.A(), indexBase, Comm);
1861 Epetra_Map overlapmap(-1, numLocNodes, (
int*)pindx.A(), indexBase, Comm);
1863 int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
1864 Gp = rcp(
new Epetra_FECrsMatrix(Copy, standardmap, 0));
1866 int* nodes =
new int[numNodesPerElem];
1867 int i=0, j=0, err=0;
1870 Epetra_SerialDenseMatrix Nodes;
1871 Epetra_SerialDenseVector Weights;
1873 int numQuadPts = Nodes.M();
1877 Epetra_SerialDenseMatrix N;
1878 N.Shape(numQuadPts,3);
1879 for (
int i=0; i<numQuadPts; i++) {
1880 N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
1881 N(i,1) = Nodes(i,0);
1882 N(i,2) = Nodes(i,1);
1886 Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
1887 Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
1889 Epetra_SerialDenseVector ly;
1890 Epetra_SerialDenseVector Nly;
1891 Epetra_SerialDenseVector lgfctn;
1892 Epetra_SerialDenseVector lgfctnNiNj;
1893 Epetra_SerialDenseMatrix lGp;
1895 ly.Size(numNodesPerElem);
1896 Nly.Size(numQuadPts);
1897 lgfctn.Size(numQuadPts);
1898 lgfctnNiNj.Size(numQuadPts);
1899 lGp.Shape(numNodesPerElem, numNodesPerElem);
1901 Epetra_SerialDenseMatrix B(2,2);
1904 for(i=0; i<numLocElems; i++) {
1906 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
1907 for (j=0; j<numNodesPerElem; j++) {
1908 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
1909 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
1913 for(
int i=0; i<2; i++) {
1914 B(i,0) = vertices(1,i)-vertices(0,i);
1915 B(i,1) = vertices(2,i)-vertices(0,i);
1920 for (j=0; j<numNodesPerElem; j++) {
1921 ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
1924 Nly.Multiply(
'N',
'N', 1.0, N, ly, 0.0);
1927 for (
int i=0; i<numNodesPerElem; i++) {
1928 for (
int j=0; j<numNodesPerElem; j++) {
1930 lGp(i,j) = adB*lgfctnNiNj.Dot(Weights);
1934 err = Gp->InsertGlobalValues(epetra_nodes, lGp, format);
1935 if (err<0)
return(err);
1940 err = Gp->GlobalAssemble();
1941 if (err<0)
return(err);
1952 void GLpApp::gpfctn(
const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv) {
1953 for (
int i=0; i<v.M(); i++) {
1954 gv(i) = 3.0*pow(v(i),2)-1.0;
1991 const Epetra_IntSerialDenseVector & ipindx,
1992 const Epetra_SerialDenseMatrix & ipcoords,
1993 const Epetra_IntSerialDenseVector & pindx,
1994 const Epetra_SerialDenseMatrix & pcoords,
1995 const Epetra_IntSerialDenseMatrix & t,
1996 const Teuchos::RCP<const Epetra_MultiVector> & y,
1997 Teuchos::RCP<Epetra_FEVector> & g)
2000 int myPID = Comm.MyPID();
2001 int numProcs = Comm.NumProc();
2003 int numLocNodes = pindx.M();
2004 int numMyLocNodes = ipindx.M();
2005 int numLocElems = t.M();
2006 int numNodesPerElem = 3;
2010 Epetra_Map standardmap(-1, numMyLocNodes, (
int*)ipindx.A(), indexBase, Comm);
2011 Epetra_Map overlapmap(-1, numLocNodes, (
int*)pindx.A(), indexBase, Comm);
2013 g = rcp(
new Epetra_FEVector(standardmap,1));
2015 int* nodes =
new int[numNodesPerElem];
2016 int i=0, j=0, err=0;
2019 Epetra_SerialDenseMatrix Nodes;
2020 Epetra_SerialDenseVector Weights;
2022 int numQuadPts = Nodes.M();
2026 Epetra_SerialDenseMatrix N;
2027 N.Shape(numQuadPts,3);
2028 for (
int i=0; i<numQuadPts; i++) {
2029 N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
2030 N(i,1) = Nodes(i,0);
2031 N(i,2) = Nodes(i,1);
2035 Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
2036 Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
2038 Epetra_SerialDenseVector ly;
2039 Epetra_SerialDenseVector Nly;
2040 Epetra_SerialDenseVector lgfctn;
2041 Epetra_SerialDenseVector lgfctnNi;
2042 Epetra_SerialDenseVector lg;
2044 ly.Size(numNodesPerElem);
2045 Nly.Size(numQuadPts);
2046 lgfctn.Size(numQuadPts);
2047 lgfctnNi.Size(numQuadPts);
2048 lg.Size(numNodesPerElem);
2050 Epetra_SerialDenseMatrix B(2,2);
2053 for(i=0; i<numLocElems; i++) {
2055 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
2056 for (j=0; j<numNodesPerElem; j++) {
2057 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
2058 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
2062 for(
int i=0; i<2; i++) {
2063 B(i,0) = vertices(1,i)-vertices(0,i);
2064 B(i,1) = vertices(2,i)-vertices(0,i);
2069 for (j=0; j<numNodesPerElem; j++) {
2070 ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
2073 Nly.Multiply(
'N',
'N', 1.0, N, ly, 0.0);
2076 for (
int i=0; i<numNodesPerElem; i++) {
2078 lg(i) = adB*lgfctnNi.Dot(Weights);
2081 err = g->SumIntoGlobalValues(epetra_nodes, lg);
2082 if (err<0)
return(err);
2087 err = g->GlobalAssemble();
2088 if (err<0)
return(err);
2100 void GLpApp::gfctn(
const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv) {
2101 for (
int i=0; i<v.M(); i++) {
2102 gv(i) = pow(v(i),3)-v(i);
2130 int MyPID = A.Comm().MyPID();
2131 int NumProc = A.Comm().NumProc();
2134 if( A.IndicesAreLocal() == false ) {
2136 std::cerr <<
"ERROR in "<< __FILE__ <<
", line " << __LINE__ << std::endl;
2137 std::cerr <<
"Function CrsMatrix2MATLAB accepts\n";
2138 std::cerr <<
"transformed matrices ONLY. Please call A.TransformToLoca()\n";
2139 std::cerr <<
"on your matrix A to that purpose.\n";
2140 std::cerr <<
"Now returning...\n";
2145 int NumMyRows = A.NumMyRows();
2150 int NumGlobalNonzeros;
2152 NumGlobalRows = A.NumGlobalRows();
2153 NumGlobalNonzeros = A.NumGlobalNonzeros();
2157 int IndexBase = A.IndexBase();
2158 if( IndexBase == 0 )
2160 else if ( IndexBase == 1)
2166 outfile <<
"A = spalloc(";
2167 outfile << NumGlobalRows <<
',' << NumGlobalRows;
2168 outfile <<
',' << NumGlobalNonzeros <<
");\n";
2171 for(
int Proc=0 ; Proc<NumProc ; ++Proc ) {
2173 if( MyPID == Proc ) {
2175 outfile <<
"\n\n% On proc " << Proc <<
": ";
2176 outfile << NumMyRows <<
" rows and ";
2177 outfile << A.NumMyNonzeros() <<
" nonzeros\n";
2180 for(
int MyRow=0 ; MyRow<NumMyRows ; ++MyRow ) {
2182 GlobalRow = A.GRID(MyRow);
2184 NumNzRow = A.NumMyEntries(MyRow);
2185 double *Values =
new double[NumNzRow];
2186 int *Indices =
new int[NumNzRow];
2188 A.ExtractMyRowCopy(MyRow, NumNzRow,
2189 NumEntries, Values, Indices);
2191 for(
int j=0 ; j<NumEntries ; ++j ) {
2192 outfile <<
"A(" << GlobalRow + IndexBase
2193 <<
"," << A.GCID(Indices[j]) + IndexBase
2194 <<
") = " << Values[j] <<
";\n";
2232 int MyPID = v.Comm().MyPID();
2233 int NumProc = v.Comm().NumProc();
2234 int MyLength = v.MyLength();
2235 int GlobalLength = v.GlobalLength();
2241 if( MyPID == 0 ) outfile <<
"v = zeros(" << GlobalLength <<
",1)\n";
2243 int NumMyElements = v.Map().NumMyElements();
2245 int * MyGlobalElements = v.Map().MyGlobalElements( );
2249 int IndexBase = v.Map().IndexBase();
2250 if( IndexBase == 0 )
2252 else if ( IndexBase == 1)
2255 for(
int Proc=0 ; Proc<NumProc ; ++Proc ) {
2257 if( MyPID == Proc ) {
2259 outfile <<
"% On proc " << Proc <<
": ";
2260 outfile << MyLength <<
" rows of ";
2261 outfile << GlobalLength <<
" elements\n";
2263 for( Row=0 ; Row<MyLength ; ++Row ) {
2264 outfile <<
"v(" << MyGlobalElements[Row] + IndexBase
2265 <<
") = " << v[Row] <<
";\n";
2300 int MyPID = v.Comm().MyPID();
2301 int NumProc = v.Comm().NumProc();
2302 int MyLength = v.MyLength();
2303 int GlobalLength = v.GlobalLength();
2309 if( MyPID == 0 ) outfile <<
"v = zeros(" << GlobalLength <<
",1)\n";
2311 int NumMyElements = v.Map().NumMyElements();
2313 int * MyGlobalElements = v.Map().MyGlobalElements( );
2317 int IndexBase = v.Map().IndexBase();
2318 if( IndexBase == 0 )
2320 else if ( IndexBase == 1)
2323 for(
int Proc=0 ; Proc<NumProc ; ++Proc ) {
2325 if( MyPID == Proc ) {
2327 outfile <<
"% On proc " << Proc <<
": ";
2328 outfile << MyLength <<
" rows of ";
2329 outfile << GlobalLength <<
" elements\n";
2331 for( Row=0 ; Row<MyLength ; ++Row ) {
2332 outfile <<
"v(" << MyGlobalElements[Row] + IndexBase
2333 <<
") = " << v[0][Row] <<
";\n";
2359 Epetra_SerialDenseMatrix & nodes,
2360 Epetra_SerialDenseVector & weights)
2373 else if (order == 2) {
2375 nodes(0,0) = (1.0-1.0/sqrt(3.0))/2.0;
2376 nodes(1,0) = (1.0+1.0/sqrt(3.0))/2.0;
2381 else if (order == 3) {
2383 nodes(0,0) = (1.0-sqrt(3.0/5.0))/2.0;
2385 nodes(2,0) = (1.0+sqrt(3.0/5.0))/2.0;
2387 weights(0) = 5.0/18.0;
2388 weights(1) = 4.0/9.0;
2389 weights(2) = 5.0/18.0;
2392 std::cout <<
"Quadrature for dim = " << dim <<
" and order = ";
2393 std::cout << order <<
" not available.\n";
2398 else if (dim == 2) {
2405 nodes(0,0) = 1.0/3.0; nodes (0,1) = 1.0/3.0;
2409 else if (order == 2) {
2411 nodes(0,0) = 1.0/6.0; nodes (0,1) = 1.0/6.0;
2412 nodes(1,0) = 2.0/3.0; nodes (1,1) = 1.0/6.0;
2413 nodes(2,0) = 1.0/6.0; nodes (2,1) = 2.0/3.0;
2415 weights(0) = 1.0/6.0;
2416 weights(1) = 1.0/6.0;
2417 weights(2) = 1.0/6.0;
2419 else if (order == 3) {
2421 nodes(0,0) = 1.0/3.0; nodes (0,1) = 1.0/3.0;
2422 nodes(1,0) = 3.0/5.0; nodes (1,1) = 1.0/5.0;
2423 nodes(2,0) = 1.0/5.0; nodes (2,1) = 3.0/5.0;
2424 nodes(3,0) = 1.0/5.0; nodes (3,1) = 1.0/5.0;
2426 weights(0) = -9.0/32.0;
2427 weights(1) = 25.0/96.0;
2428 weights(2) = 25.0/96.0;
2429 weights(3) = 25.0/96.0;
2432 std::cout <<
"Quadrature for dim = " << dim <<
" and order = ";
2433 std::cout << order <<
" not available.\n";
2439 std::cout <<
"Quadrature for dim = " << dim <<
" not available.\n";
Teuchos::RCP< Epetra_FECrsMatrix > getR()
Teuchos::RCP< Epetra_FEVector > getq()
Teuchos::RCP< const Epetra_IntSerialDenseVector > getpindx()
int assemblyFECrs(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, Teuchos::RCP< Epetra_FECrsMatrix > &, Teuchos::RCP< Epetra_FEVector > &, bool)
void PrintVec(const Teuchos::RCP< const Epetra_Vector > &x)
Outputs the solution vector to files.
Epetra_SerialDenseMatrix NNw
Teuchos::RCP< const Epetra_IntSerialDenseMatrix > gett()
void computeNpy(const Teuchos::RCP< const Epetra_MultiVector > &y)
Calls the function that computes the Jacobian of the nonlinear term.
Epetra_SerialDenseMatrix * NdNdx2Nw
int nonlinjac(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Teuchos::RCP< const Epetra_MultiVector > &, Teuchos::RCP< Epetra_FECrsMatrix > &)
Epetra_SerialDenseMatrix S3
void computeNy(const Teuchos::RCP< const Epetra_MultiVector > &y)
Calls the function that computes the nonlinear term.
void gfctn(const Epetra_SerialDenseVector &, Epetra_SerialDenseVector &)
bool Vector2MATLAB(const Epetra_Vector &, std::ostream &)
bool FEVector2MATLAB(const Epetra_FEVector &, std::ostream &)
int nonlinvec(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Teuchos::RCP< const Epetra_MultiVector > &, Teuchos::RCP< Epetra_FEVector > &)
Epetra_SerialDenseMatrix S1
Teuchos::RCP< const Epetra_IntSerialDenseMatrix > gete()
int compproduct(Epetra_SerialDenseVector &, double *, double *, double *)
int meshreader(const Epetra_Comm &, Epetra_IntSerialDenseVector &, Epetra_SerialDenseMatrix &, Epetra_IntSerialDenseVector &, Epetra_SerialDenseMatrix &, Epetra_IntSerialDenseMatrix &, Epetra_IntSerialDenseMatrix &, const char geomFileBase[], const bool trace=true, const bool dumpAll=false)
double determinant(const Epetra_SerialDenseMatrix &)
Provides the interface to generic abstract vector libraries.
void Print(std::ostream &os) const
Teuchos::RCP< Epetra_FECrsMatrix > getNpy()
Transform to form the explicit transpose of a Epetra_RowMatrix.
void gpfctn(const Epetra_SerialDenseVector &v, Epetra_SerialDenseVector &gv)
Epetra_SerialDenseMatrix S2
void rect2DMeshGenerator(const int numProc, const int procRank, const double len_x, const double len_y, const int local_nx, const int local_ny, const int bndy_marker, Epetra_IntSerialDenseVector *ipindx_out, Epetra_SerialDenseMatrix *ipcoords_out, Epetra_IntSerialDenseVector *pindx_out, Epetra_SerialDenseMatrix *pcoords_out, Epetra_IntSerialDenseMatrix *t_out, Epetra_IntSerialDenseMatrix *e_out, std::ostream *out, const bool dumpAll)
Generate a simple rectangular 2D triangular mesh that is only partitioned between processors in the y...
Epetra_SerialDenseMatrix * NdNdx1Nw
void computeAugmat()
Assembles the augmented system (KKT-type) matrix.
Teuchos::RCP< const Epetra_SerialDenseMatrix > getpcoords()
int assemble_bdry(const Epetra_Comm &Comm, const Epetra_IntSerialDenseVector &ipindx, const Epetra_SerialDenseMatrix &ipcoords, const Epetra_IntSerialDenseVector &pindx, const Epetra_SerialDenseMatrix &pcoords, const Epetra_IntSerialDenseMatrix &t, const Epetra_IntSerialDenseMatrix &e, Teuchos::RCP< Epetra_FECrsMatrix > *B, Teuchos::RCP< Epetra_FECrsMatrix > *R)
void g2pfctn(const Epetra_SerialDenseVector &, Epetra_SerialDenseVector &)
int compproduct(Epetra_SerialDenseVector &, double *, double *)
int lassembly(const Epetra_SerialDenseMatrix &, const Epetra_SerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_SerialDenseVector &, const Epetra_SerialDenseVector &, const Usr_Par &, Epetra_SerialDenseMatrix &, Epetra_SerialDenseVector &)
Teuchos::RCP< const Epetra_IntSerialDenseVector > getipindx()
Teuchos::RCP< Epetra_FECrsMatrix > getH()
Teuchos::RCP< const Epetra_SerialDenseMatrix > getipcoords()
void computeAll(const GenSQP::Vector &x)
Calls functions to compute nonlinear quantities and the augmented system matrix.
Teuchos::RCP< Epetra_FECrsMatrix > getB()
int nonlinhessvec(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Teuchos::RCP< const Epetra_MultiVector > &, const Teuchos::RCP< const Epetra_MultiVector > &, const Teuchos::RCP< const Epetra_MultiVector > &, Teuchos::RCP< Epetra_FEVector > &)
Teuchos::RCP< Epetra_FECrsMatrix > getA()
Teuchos::RCP< Epetra_CrsMatrix > getAugmat()
bool CrsMatrix2MATLAB(const Epetra_CrsMatrix &, std::ostream &)
Epetra_SerialDenseVector Nw
Epetra_SerialDenseMatrix * NNNw
int solveAugsys(const Teuchos::RCP< const Epetra_MultiVector > &rhsy, const Teuchos::RCP< const Epetra_MultiVector > &rhsu, const Teuchos::RCP< const Epetra_MultiVector > &rhsp, const Teuchos::RCP< Epetra_MultiVector > &y, const Teuchos::RCP< Epetra_MultiVector > &u, const Teuchos::RCP< Epetra_MultiVector > &p, double tol)
Solves augmented system.
Teuchos::RCP< Epetra_FEVector > getNy()
int quadrature(const int, const int, Epetra_SerialDenseMatrix &, Epetra_SerialDenseVector &)
std::ostream & operator<<(std::ostream &, const Usr_Par &)
The GenSQP::Vector / (y,u) Epetra_MultiVector adapter class.
Teuchos::RCP< const Epetra_Comm > getCommPtr()
Teuchos::RCP< Epetra_FEVector > getb()
int assemble(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, Teuchos::RCP< Epetra_FECrsMatrix > &, Teuchos::RCP< Epetra_FECrsMatrix > &, Teuchos::RCP< Epetra_FEVector > &)
int inverse(const Epetra_SerialDenseMatrix &, Epetra_SerialDenseMatrix &)