Index: ../trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp (revision 13938) +++ ../trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp (revision 13939) @@ -1,47 +0,0 @@ -/*!\file CreateJacobianMatrixx - * \brief: create system matrices (stiffness matrix, loads vector) - */ - -#include "./CreateJacobianMatrixx.h" -#include "../../shared/shared.h" -#include "../../include/include.h" -#include "../../io/io.h" -#include "../../toolkits/toolkits.h" -#include "../../EnumDefinitions/EnumDefinitions.h" - -void CreateJacobianMatrixx(Matrix** pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,IssmDouble kmax){ - - int i,connectivity; - int numberofdofspernode; - int fsize,configuration_type; - Element *element = NULL; - Load *load = NULL; - Matrix* Jff = NULL; - - /*Checks*/ - _assert_(nodes && elements); - - /*Recover some parameters*/ - parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum); - fsize=nodes->NumberOfDofs(configuration_type,FsetEnum); - numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum); - - /*Initialize Jacobian Matrix*/ - Jff=new Matrix(fsize,fsize,connectivity,numberofdofspernode); - - /*Create and assemble matrix*/ - for(i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - element->CreateJacobianMatrix(Jff); - } - for (i=0;iSize();i++){ - load=(Load*)loads->GetObjectByOffset(i); - if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff); - if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax); - } - Jff->Assemble(); - - /*Assign output pointer*/ - *pJff=Jff; -} Index: ../trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h =================================================================== --- ../trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h (revision 13938) +++ ../trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h (revision 13939) @@ -1,14 +0,0 @@ -/*!\file: CreateJacobianMatrixx.h - * \brief header file for degree of freedoms distribution routines. - */ - -#ifndef _CREATEJACOBIANMATRIXX_H -#define _CREATEJACOBIANMATRIXX_H - -#include "../../Container/Container.h" -#include "../../classes/objects/objects.h" - -/* local prototypes: */ -void CreateJacobianMatrixx(Matrix** pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,IssmDouble kmax); - -#endif /* _CREATEJACOBIANMATRIXX_H */ Index: ../trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 13938) +++ ../trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 13939) @@ -14,7 +14,7 @@ int configuration_type; int fsize; IssmDouble *ug_serial = NULL; - bool oldalloc = false; + bool oldalloc = true; /*first figure out fsize: */ parameters->FindParam(&configuration_type,ConfigurationTypeEnum); Index: ../trunk-jpl/src/c/modules/modules.h =================================================================== --- ../trunk-jpl/src/c/modules/modules.h (revision 13938) +++ ../trunk-jpl/src/c/modules/modules.h (revision 13939) @@ -90,7 +90,6 @@ #include "./Solverx/Solverx.h" #include "./SpcNodesx/SpcNodesx.h" #include "./SurfaceAreax/SurfaceAreax.h" -#include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h" #include "./TriaSearchx/TriaSearchx.h" #include "./TriMeshx/TriMeshx.h" #include "./TriMeshProcessRiftsx/TriMeshProcessRiftsx.h" Index: ../trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.cpp (revision 13938) +++ ../trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.cpp (revision 13939) @@ -19,6 +19,7 @@ int Kfsm,Kfsn; int global_m,global_n; bool fromlocalsize = true; + bool oldalloc = true; if(VerboseModule()) _pprintLine_(" Dirichlet lifting applied to load vector"); @@ -30,7 +31,11 @@ /*pf = pf - Kfs * y_s;*/ Kfs->GetLocalSize(&Kfsm,&Kfsn); - Kfsy_s=new Vector(Kfsm,fromlocalsize); + if(oldalloc) + Kfsy_s=new Vector(Kfsm,fromlocalsize); + else + Kfsy_s=new Vector(Kfsm,global_m); + if (flag_ys0){ /*Create y_s0, full of 0: */ Index: ../trunk-jpl/src/c/Makefile.am =================================================================== --- ../trunk-jpl/src/c/Makefile.am (revision 13938) +++ ../trunk-jpl/src/c/Makefile.am (revision 13939) @@ -298,8 +298,6 @@ ./modules/EnumToStringx/EnumToStringx.h\ ./modules/StringToEnumx/StringToEnumx.cpp\ ./modules/StringToEnumx/StringToEnumx.h\ - ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp\ - ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h\ ./modules/ConstraintsStatex/ConstraintsStatex.cpp\ ./modules/ConstraintsStatex/ConstraintsStatex.h\ ./modules/ConstraintsStatex/ConstraintsStateLocal.h\ Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 13938) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 13939) @@ -352,7 +352,7 @@ Vector *pf = NULL; Vector *df = NULL; - bool oldalloc=false; + bool oldalloc=true; /*retrieve parameters: */ this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); @@ -537,9 +537,9 @@ xDelete(next_l); /*sum over all cpus*/ - #ifdef _HAVE_MPI_ +#ifdef _HAVE_MPI_ MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,MPI_INT,MPI_SUM,IssmComm::GetComm()); - #endif +#endif xDelete(connectivity_clone); if(set1enum==FsetEnum){ @@ -573,6 +573,44 @@ *po_nnz=o_nnz; }/*}}}*/ +void FemModel::CreateJacobianMatrixx(Matrix** pJff,IssmDouble kmax){/*{{{*/ + + int i,connectivity; + int numberofdofspernode; + int fsize,configuration_type; + Element *element = NULL; + Load *load = NULL; + Matrix* Jff = NULL; + + /*Checks*/ + _assert_(nodes && elements); + + /*Recover some parameters*/ + parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum); + fsize=nodes->NumberOfDofs(configuration_type,FsetEnum); + numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum); + + /*Initialize Jacobian Matrix*/ + this->AllocateSystemMatrices(&Jff,NULL,NULL,NULL); + Jff=new Matrix(fsize,fsize,connectivity,numberofdofspernode); + + /*Create and assemble matrix*/ + for(i=0;iSize();i++){ + element=dynamic_cast(elements->GetObjectByOffset(i)); + element->CreateJacobianMatrix(Jff); + } + for (i=0;iSize();i++){ + load=(Load*)loads->GetObjectByOffset(i); + if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff); + if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax); + } + Jff->Assemble(); + + /*Assign output pointer*/ + *pJff=Jff; + +}/*}}}*/ int FemModel::UpdateVertexPositionsx(void){ /*{{{*/ int i; Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 13938) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 13939) @@ -52,6 +52,7 @@ /*Methods:*/ void AllocateSystemMatrices(Matrix** pKff,Matrix** pKfs,Vector** pdf,Vector** ppf); + void CreateJacobianMatrixx(Matrix** pJff,IssmDouble kmax); void Echo(); void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels); void MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum);