2828
2929
3030#include < Core/Algorithms/Legacy/FiniteElements/BuildRHS/BuildFEVolRHS.h>
31+ #include < Core/Algorithms/Base/AlgorithmPreconditions.h>
32+ #include < Core/Algorithms/Base/AlgorithmVariableNames.h>
3133#include < Core/Datatypes/DenseMatrix.h>
3234#include < Core/Datatypes/Legacy/Field/Mesh.h>
3335#include < Core/Datatypes/Legacy/Field/VMesh.h>
@@ -48,25 +50,223 @@ using namespace SCIRun::Core::Algorithms;
4850using namespace SCIRun ::Core::Algorithms::FiniteElements;
4951// using namespace SCIRun::Core::Logging;
5052
53+ class FEMVolRHSBuilder
54+ {
55+ public:
56+
57+ FEMVolRHSBuilder (const AlgorithmBase *algo) :
58+ algo_ (algo),
59+ barrier_ (" FEMVolRHSBuilder Barrier" , numprocessors_)
60+ {
61+ }
62+
63+ bool build_RHS_Vol (FieldHandle input, MatrixHandle vtable, MatrixHandle& output);
64+
65+ private:
66+ const AlgorithmBase *algo_;
67+ Barrier barrier_;
68+
69+ VMesh* mesh_;
70+ VField *field_;
71+
72+ MatrixHandle rhsmatrixhandle_;
73+ DenseMatrix* rhsmatrix_;
74+
75+ // num processors
76+ int numprocessors_;
77+
78+ std::vector<bool > success_;
79+
80+ index_type domain_dimension;
81+
82+ index_type local_dimension_nodes;
83+ index_type local_dimension_add_nodes;
84+ index_type local_dimension_derivatives;
85+ index_type local_dimension;
86+
87+ index_type global_dimension_nodes;
88+ index_type global_dimension_add_nodes;
89+ index_type global_dimension_derivatives;
90+ index_type global_dimension;
91+
92+ bool use_vector_;
93+ std::vector<std::pair<std::string, Vector> > vectors_;
94+
95+ bool use_scalars_;
96+ std::vector<std::pair<std::string, double > > scalars_;
97+
98+ // Entry point for the parallel version
99+ void parallel (int proc);
100+
101+ private:
102+
103+ void create_numerical_integration (std::vector<VMesh::coords_type > &p,
104+ std::vector<double > &w,
105+ std::vector<std::vector<double > > &d);
106+ bool build_local_matrix (VMesh::Elem::index_type c_ind,
107+ index_type row,
108+ double &l_val,
109+ std::vector<VMesh::coords_type> &p,
110+ std::vector<double > &w,
111+ std::vector<std::vector<double > > &d);
112+ bool build_local_matrix_regular (VMesh::Elem::index_type c_ind,
113+ index_type row,
114+ double &l_val,
115+ std::vector<VMesh::coords_type> &p,
116+ std::vector<double > &w,
117+ std::vector<std::vector<double > > &d,
118+ std::vector<std::vector<double > > &precompute);
119+ bool setup ();
120+
121+ };
122+
51123AlgorithmInputName BuildFEVolRHSAlgo::Mesh (" Mesh" );
52124AlgorithmInputName BuildFEVolRHSAlgo::Vector_Table (" Vector_Table" );
53125AlgorithmOutputName BuildFEVolRHSAlgo::RHS (" RHS" );
126+ AlgorithmParameterName BuildFEVolRHSAlgo::vectorTableBasisMatrices () { return AlgorithmParameterName (" vectorTableBasisMatrices" ); }
54127
55128BuildFEVolRHSAlgo::BuildFEVolRHSAlgo ()
56129{
57-
130+ addParameter ( vectorTableBasisMatrices (), false );
58131}
59132
60-
61133DenseMatrixHandle BuildFEVolRHSAlgo::run (FieldHandle input, DenseMatrixHandle ctable) const
62134{
63- DenseMatrixHandle a;
64-
65-
135+ DenseMatrixHandle rhs_matrix;
66136
137+ if (input)
138+ {
139+ THROW_ALGORITHM_INPUT_ERROR (" Could not obtain input field" );
140+ }
141+
142+ if (!input->vfield ()->is_vector ())
143+ {
144+ THROW_ALGORITHM_INPUT_ERROR (" This function is only defined for elements with vector data" );
145+ }
146+
147+ if (input->vfield ()->basis_order ()!=0 )
148+ {
149+ THROW_ALGORITHM_INPUT_ERROR (" This function has only been defined for data that is located at the elements" );
150+ }
151+
152+ if (ctable)
153+ {
154+ if ((ctable->ncols () != 1 )&&(ctable->ncols () != 3 ))
155+ {
156+ THROW_ALGORITHM_INPUT_ERROR (" Vector table needs to have 1 or 3" );
157+
158+ }
159+ if (ctable->nrows () == 0 )
160+ {
161+ THROW_ALGORITHM_INPUT_ERROR (" Vector table is empty" );
162+ }
163+ }
164+
165+ FEMVolRHSBuilder builder (this );
67166
167+ if (get (vectorTableBasisMatrices ()).toBool ())
168+ {
169+ #ifdef SCIRUN4_CODE_TO_BE_ENABLED_LATER
170+ BuildFEVolRHSPrivateData* privatedata;
171+ get_privatedata (privatedata);
172+ #endif
173+
174+ if (ctable)
175+ {
176+ std::vector<std::pair<std::string,Tensor> > tens;
177+ #ifdef SCIRUN4_CODE_TO_BE_ENABLED_LATER
178+ input->get_property (" conductivity_table" ,tens);
179+ if (tens.size () > 0 )
180+ {
181+ vtable = new DenseMatrix (tens.size (),1 );
182+ double * data = vtable->get_data_pointer ();
183+ for (size_t i=0 ; i<tens.size ();i++)
184+ {
185+ double t = tens[i].second .mat_ [0 ][0 ];
186+ data[i] = t;
187+ }
188+ }
189+ #endif
190+ }
191+
192+ if (ctable)
193+ {
194+ size_type nconds = ctable->nrows ();
195+ if ((input->vmesh ()->generation () != generation_)||
196+ (basis_fevolrhs_))
197+ {
198+ /* MatrixHandle con = new DenseMatrix(nconds,1);
199+ double* data = con->get_data_pointer();
200+ for (size_type i=0; i<nconds;i++) data[i] = 0.0;
201+ builder->build_RHS_Vol(input,con,privatedata->basis_fevolrhs_);
202+ if (privatedata->basis_fevolrhs_.get_rep() == 0)
203+ {
204+ error("Failed to build FEVolRHS structure");
205+ algo_end(); return(false);
206+ }
207+ privatedata->basis_values_.resize(nconds);
208+ for (size_type s=0; s< nconds; s++)
209+ {
210+ MatrixHandle temp;
211+ data[s] = 1.0;
212+ builder->build_RHS_Vol(input,con,temp);
213+ if (temp.get_rep() == 0)
214+ {
215+ error("Failed to build FE component for one of the tissue types");
216+ algo_end(); return(false);
217+ }
218+
219+ SparseRowMatrix *m = temp->sparse();
220+ privatedata->basis_values_[s].resize(m->get_nnz());
221+ for (size_type p=0; p< m->get_nnz(); p++)
222+ {
223+ privatedata->basis_values_[s][p] = m->get_value(p);
224+ }
225+ data[s] = 0.0;
226+ }
227+
228+ privatedata->generation_ = input->vmesh()->generation();
229+ */
230+ }
231+
232+ /*
233+ output = privatedata->basis_fevolrhs_;
234+ output.detach();
235+
236+ SparseRowMatrix *m = output->sparse();
237+ double *sum = m->get_vals();
238+ double *cdata = vtable->get_data_pointer();
239+ size_type n = vtable->ncols();
240+
241+ if (privatedata->basis_values_.size() > 0)
242+ for (size_t p=0; p < privatedata->basis_values_[0].size(); p++) sum[p] = 0.0;
243+
244+ for (int s=0; s<nconds; s++)
245+ {
246+ double weight = cdata[s*n];
247+ for (size_t p=0; p < privatedata->basis_values_[s].size(); p++)
248+ {
249+ sum[p] += weight * privatedata->basis_values_[s][p];
250+ }
251+ }
252+ */
253+ }
254+ else
255+ {
256+ THROW_ALGORITHM_INPUT_ERROR (" No vector table present: The generate_basis option only works for indexed vectors" );
257+ }
258+
259+ }
260+
261+ // builder->build_RHS_Vol(input,vtable,output);
262+ /*
263+ if (output)
264+ {
265+ THROW_ALGORITHM_INPUT_ERROR("Could not build output matrix");
266+ }
267+ */
68268
69- return a ;
269+ return rhs_matrix ;
70270}
71271
72272AlgorithmOutput BuildFEVolRHSAlgo::run_generic (const AlgorithmInput& input) const
0 commit comments